Permutation inference for canonical correlation analysis.

Canonical correlation analysis (CCA) has become a key tool for population neuroimaging, allowing investigation of associations between many imaging and non-imaging measurements. As age, sex and other variables are often a source of variability not of direct interest, previous work has used CCA on residuals from a model that removes these effects, then proceeded directly to permutation inference. We show that a simple permutation test, as typically used to identify significant modes of shared variation on such data adjusted for nuisance variables, produces inflated error rates. The reason is that residualisation introduces dependencies among the observations that violate the exchangeability assumption. Even in the absence of nuisance variables, however, a simple permutation test for CCA also leads to excess error rates for all canonical correlations other than the first. The reason is that a simple permutation scheme does not ignore the variability already explained by previous canonical variables. Here we propose solutions for both problems: in the case of nuisance variables, we show that transforming the residuals to a lower dimensional basis where exchangeability holds results in a valid permutation test; for more general cases, with or without nuisance variables, we propose estimating the canonical correlations in a stepwise manner, removing at each iteration the variance already explained, while dealing with different number of variables in both sides. We also discuss how to address the multiplicity of tests, proposing an admissible test that is not conservative, and provide a complete algorithm for permutation inference for CCA. Published by Elsevier Inc.

Introduction

Canonical correlation analysis (cca) (Jordan, 1875; Hotelling, 1936) is a multivariate method that aims at reducing the correlation structure between two sets of variables to the simplest possible form (hence the name ‘‘canonical’’) through linear transformations of the variables within each set. Put simply, given two sets of variables, the method seeks linear mixtures within each set, such that each resulting mixture from one set is maximally correlated with a corresponding mixture from the other set, but uncorrelated with all other mixtures in either set.

From a peak use through from the late 1970’s until mid-1980’s, the method has recently regained popularity, presumably thanks to its ability to uncover latent, common factors underlying association between multiple measurements obtained, something relevant in recent research, particularly in the field of brain imaging, that uses high dimensional phenotyping and investigates between-subject variability across multiple domains. This is in contrast to initial studies that introduced cca to the field (Friston et al., 1995, 1996; Worsley, 1997; Friman et al., 2001, 2002, 2003) for investigation of signal variation in functional magnetic resonance imaging (fmri) time series. For example, Smith et al. (2015) used cca to identify underlying factors associating brain connectivity features to various demographic, psychometric, and lifestyle measures; Rosa et al. (2015) used a sparse cca method to investigate differences in brain perfusion after administration of two distinct antipsychotic drugs; Miller et al. (2016) used cca to identify associations between imaging and non-imaging variables in the uk Biobank; Drysdale et al. (2017) used cca to investigate associations between brain connectivity and clinical assessments, and found two canonical variables that would allow classification of participants into distinct categories (but see Dinga et al., 2019); Kernbach et al. (2018) used cca to identify connectivity patterns in the default mode network associated with patterns of connectivity elsewhere in the brain; Bijsterbosch et al. (2018), Xia et al. (2018), and Mihalik et al. (2019) likewise used cca to identify associations between functional connectivity and various indices of behaviour and psychopathology, whereas Sui et al. (2018) used a combination of multivariate methods, including cca, to investigate brain networks associated with composite cognitive scores; Li et al. (2019) used cca to investigate, among subjects, the topography of the global fmri signal and its relationship with a number of cognitive and behavioral measurements; Ing et al. (2019) used cca to identify symptom groups that were correlated with brain regions assessed through a diverse set of imaging modalities; Alnæs et al. (2020) used cca to investigate the association between imaging measurements and cognitive, behavioral, psychosocial and socioeconomic indices; Clemens et al. (2020) used a combination of pattern classification algorithms and cca to study imaging and behavioral correlates of the subjective perception of oneself belonging to a particular gender. In most of these between-subject, group level studies, putative nuisance variables or confounds were regressed out from the data before proceeding to inference, and all of them used some form of permutation test to assess the significance of the results. In a recent review, Wang et al. (2020) described a permutation procedure for cca as ‘‘a random shuffling of the rows (…) of the two variable sets’’.

Permutation tests are well known and widely used. Among their many benefits, these tests lead to valid inferences while requiring assumptions that are commonly satisfied in between-subject analyses, such that of exchangeability of observations under the null hypothesis. However, here we show that simple implementations of permutation inference for cca are inadequate on four different grounds. First, simple, uncorrected permutation p-values are not guaranteed to be monotonically related to the canonical correlations, leading to inadmissible results; for the same reason, multiple testing correction using false discovery rate is also inadmissible. Second, except for the highest canonical correlation, a simple one-step estimation of all others without considering the variability already explained by previous canonical variables in relation to each of them also leads to inflated per comparison error rates and thus, invalid results. Third, regressing out nuisance variables without consideration to the introduction of dependencies among observations caused by residualisation leads to an invalid test, with excess false positives. Fourth, multiple testing correction using the distribution of the maximum test statistic leads to conservative results, except for the highest canonical correlation.

In this paper we explain and discuss in detail each of these problems, and offer solutions that address each of them. In particular, we propose a stepwise estimation method for the canonical correlations and canonical variables that remains valid even when the number of variables is not the same for both sides of cca. We propose a method that transforms residualised data to a lower dimensional basis where exchangeability — as required for the validity of permutation tests — holds. We also propose that inference that considers multiple canonical correlations should use a closed testing procedure that is more powerful than the usual correction method used in permutation tests that use the distribution of the maximum statistic; the procedure also ensures a monotonic relationship between p-values and canonical correlations. Finally, we provide a complete, general algorithm for valid inferences for cca.

Theory

Notation and general aspects

Thorough definition and derivation of cca can be found in many classical textbooks on multivariate analysis (e.g., Kendall, 1975; Mardia et al., 1979; Brillinger, 1981; Muirhead, 1982; Seber, 1984; Krzanowski, 1988; Anderson, 2003); the reader is referred to these for a comprehensive overview. Here we present concisely and only allude to the distinction between population (ρ) and sample (r) canonical correlations where strictly needed. Let and be each one a collection of, respectively, P and Q variables observed from N subjects, . Without loss of generality, assume that , that the columns of and are mean-centered, that these matrices are of full rank, and define , , and . The goal of cca is to identify canonical coefficients or canonical weights and , , such that the pairs of canonical variables, defined as:have correlations that are maximal, under the constraint that . Estimation of and amounts to finding the K solutions to:where the unknowns are , and ; are the sample canonical correlations, i.e., the correlations between the estimated canonical variables and . The coefficients are eigenvectors of , whereas are eigenvectors of ; the respective eigenvalues (for either or , as these eigenvalues are the same) are . For convenience, we call canonical component the ensemble formed by the k-th canonical correlation, its corresponding pair of canonical variables, and associated pair of canonical coefficients; canonical variables may also be termed modes of variation.

The typical method for estimation involves an iterative procedure that finds one and at a time, with computed as a function of these. However, the method proposed by Bjorck and Golub (1973) is more concise and numerically more stable; it is described in the Appendix (Algorithm 3). The canonical correlations are then produced in descending order, ; this positiveness of all canonical correlations is a consequence of these values being explicitly maximised during estimation; reversal of the sign of the coefficients can always be accompanied by the reversal of the sign of the corresponding coefficients in the other side (and of and ), to no net effect on . That is, the signs of the canonical variables and coefficients are indeterminate, and any solution is arbitrary; nothing can be concluded about the specific direction of effects with cca.

Parametric inference

The distribution of the sample canonical correlations is intractable, even under assumptions of normality and independence among subjects, and is a function of the population correlations (Constantine, 1963; James, 1964). This difficulty led to the development of a rich asymptotic theory (Fisher, 1939; Hsu, 1941; Lawley, 1959; Fujikoshi, 1977; Glynn and Muirhead, 1978). However, these approximations have been shown to be extremely sensitive to departures from normality, or require additional terms that further complicate their use (Muirhead and Waternaux, 1980); Brillinger (1981) recommended resampling methods to estimate parameters used by normal approximations, which otherwise can be biased (Anderson, 2003). These difficulties pose challenges for inference. Even though some computationally efficient algorithms have been proposed (Koev and Edelman, 2006), these tests continue to be rarely used.

Instead, a test based on whether a certain number of correlations are equal to zero has been proposed. The null hypothesis is , i.e., , for , that is, the null is that population canonical correlations (the smaller ones) are zero (Bartlett, 1938, 1947; Marriott, 1952; Lawley, 1959; Fujikoshi, 1974), versus the alternative that at least one is not, i.e., . The test is based on the statistic proposed by Wilks (1935), as:where the constant if there are no nuisance variables (Section 2.6). Under the null hypothesis, follows an approximate distribution with degrees of freedom if each of the columns of and have values that are independent and identically distributed following a normal distribution (but see Glynn and Muirhead, 1978, for a different expression). Unfortunately, this test is sensitive to departures from normality, particularly in the presence of outliers, and its use has been questioned (Seber, 1984; Harris, 1976).

Another test statistic is based on Roy (1953)1 and is simply:

The critical values for the corresponding parametric distribution at a given test level α can be found in the charts provided by Heck (1960), using as parameters , , and (Lee, 1978), where the constant if there are no nuisance variables (Section 2.6), or in tables provided by Kres (1975); more recent approximations for normally distributed data can be found in Chiani (2016) and Johnstone and Nadler (2017). Some approximations, however, are considered conservative (Harris, 1976, 2013). Note that, while Roy (1953) proposed the use of the largest value as test statistic, which would then be , any given null at position k must have already considered the previous canonical components, from 1 until , such that the maximum statistic, after the previous canonical correlations have been removed from the model, is always the current one. A similar reasoning holds for the smallest canonical correlations in the test proposed by Wilks (1935). This feature is exploited in the stepwise rejective procedure proposed in Section 2.5.

Permutation inference

The above problems can be eschewed with the use of resampling methods, such as permutation. An intuitive (but inadequate) permutation test for cca could be constructed by randomly permuting the rows of or . For each shuffling of the data, indicated by , a new set of canonical correlations and associated statistics would be computed. A p-value would be obtained as , where is the indicator (Kronecker) function, which evaluates as 1 if the condition inside the brackets is true, or 0 otherwise, and the index corresponds to the unpermuted data (i.e., no permutation, with the data in their original ordering).

Such a naïve procedure, however, would ignore the fact that, this resampling scheme matches the first null hypothesis , but not the subsequent ones. For a given canonical correlation at position k being tested, , one must generate a permutation that satisfies the corresponding null , but not necessarily . Otherwise, latent effects represented by the corresponding earlier canonical variables and would, in the procedure above, remain in the and at the time these are permuted. However, the variance associated with these earlier canonical variables would have already been explained through the rejection of their respective null hypothesis up to . This variance is now a nuisance for the positions from k (inclusive) onward. It contains information that are not pertinent to position k and subsequent ones, and that therefore should not be used to build the null distribution. That is variance should not be re-used in the shufflings that lead to a given or subsequent correlations.

Fortunately, cca is invariant to linear transformations that mix the variables in or in . Since the canonical variables are themselves linear transformations of these input variables (Equation (1)), a second cca using and in place of the initial and leads to the same solutions. Yet, unless , will not span the same space as . In principle, this would be inconsequential as far as the canonical variables are concerned. However, ignoring the variability in not contained in would again affect the p-values should and be used in a permutation test, as the permuted data would not be representative of the original (unpermuted) that led to these initial canonical variables. To mitigate the problem, include into the matrix of canonical coefficients their orthogonal complement, i.e., compute , then use instead of as a replacement for . In this paper we adopted the convention that , but the same procedure works in reverse and, algorithmically, it might as well be simpler to compute also a and use it instead of as a replacement for . If , then .

While these transformations do not change in any way the canonical components, they allow the construction of an improved algorithm that addresses the issue of variability already explained by canonical variables of lower rank (i.e., the ones with order indices smaller than that of a given one). It consists of running an initial cca using and to obtain and , then subject these to a second cca and permutation testing while, crucially, at each permutation, iteratively repeating cca K times, each using not the whole and , but only from the k-th component onwards, i.e., and for the test about the k-th canonical correlation. Of note, a test level α does not need to be specified at the time in which the above iterative (stepwise) procedure is performed; instead, and in combination with the multiple testing procedure described below, adjusted p-values are are computed, which then are used to accept or reject the null for the k-th correlation. Algorithm 1 (Section 2.8) shows the procedure in detail (the algorithm contains other details that are discussed in the next sections).

A number of further aspects need be considered in permutation tests: the number of possible reorderings of the data, the need for permutations that break the association between the variables being tested, the random selection of permutations from the permutation set when not all possible permutations can be used, the choice of the test statistic, how to correct for the multiplicity of tests, the number of permutations to allow narrow confidence intervals around p-values, among others. These topics have been discussed in Winkler et al. (2014, 2016) and references therein and will not be repeated here. However, for cca, some aspects deserve special treatment and are considered below.

Choice of the statistic

Asymptotically, using Wilks’ statistic or Roy’s are expected to lead to the same conclusion since all correlations are sorted in descending order: if , then all subsequent ones must be zero; likewise, if , then clearly at least one correlation between k and K is larger than zero, which has to include itself. Moreover, permutation under the null is justifiable in the complete absence of association between the two sets, which implies that, under the null , all correlations are equal to zero. With finite data, however, one statistic can be more powerful than the other in different settings; their relative performance is studied in Sections 3, 4.

Computationally, Wilks’ requires more operations to be performed compared to Roy’s statistic. Since the relationship between and is monotonic, the two are permutationally equivalent, and using alone is sufficient, which makes Roy’s the absolute fastest. However, even for Wilks’, the amount of computation required is negligible compared to the overall number of operations needed for estimation of the canonical coefficients, such that in practice, the choice between the two should be in terms of power.

In either case, while inference refers to the respective null hypothesis at position k, it is not to be understood as inference on the index k. Rather, the null is merely conditional on the nulls for all previous correlations from 1 to having been rejected. Rejecting the null implies that the correlation observed at position k is too high under the null hypothesis of no association between the two variable sets after all previous (from 1 to ) canonical variables have been sequentially removed, as described in Section 2.3. In Algorithm 1 (Section 2.8) this is done in line 29, that uses as inputs to cca the precomputed canonical variables only from position k onwards, as opposed to all of them.

Multiplicity of tests

For either of these two test statistics, the ordering of the canonical correlations from larger (farther from zero) to smaller (closer to zero) imply that rejection of the null hypothesis at each k must happen sequentially, starting from , using the respective test statistic and associated p-value until the null , for some , is not rejected at a predefined test level α. Then, at that position k, the procedure stops, and the null is retained from that position (inclusive) onward until the final index K.

The ordering of the canonical correlations brings additional consequences. First, because rejection of implies rejection of all joint (intersection) hypotheses that include , that is , such sequentially rejective procedure is also a closed testing procedure (ctp), which controls the amount of any type i error across all tests, i.e., the familywise error rate (fwer) in the strong sense (Marcus et al., 1976; Hochberg and Tamhane, 1987). Thus, there is no need for further correction for multiple testing. Another way of stating the same is that the test for a given , , has been ‘‘protected’’ by the test at the position at the level α. Adjusted p-values (in the fwer sense) can be computed as , that is, the fwer-adjusted p-value for the canonical component k is the cumulative maximum p-value up to position k. Such adjusted p-values can be considered significant if their value is below the desired test level α.

The second consequence is that fwer adjustment of p-values using the distribution of the maximum statistic (not to be confused with the cumulative maximum described in the above paragraph) will be conservative for all canonical components except the first. The reason is that the maximum statistic is always the distribution of the first, which is stochastically dominant over all others. The distribution of the maximum is usually sought as a shortcut to a ctp when the condition of subset pivotality holds (Westfall and Young, 1993), as that reduces the computational burden from tests to only K tests. Interestingly, the ordering of the canonical correlations from largest to smallest leads to a ctp that does not use the distribution of the maximum, and yet requires only K tests, regardless of whether subset pivotality holds.

A third consequence is that using permutation p-values outside the above sequentially rejective procedure that controls fwer is not appropriate since these simple, uncorrected p-values are not guaranteed to be monotonically related to the canonical correlations . Attempts to use these uncorrected p-values outside a ctp would lead to paradoxical results whereby higher, stronger canonical correlations might not be considered significant, yet later ones, smaller, weaker, could be so; that is, it would make the test inadmissible (Lehmann and Romano, 2005, p. 232). For the same reason, such simple p-values should not be subjected to correction using false discovery rate (fdr; Benjamini and Hochberg, 1995), because the ordering of p-values for fdr, from smallest to largest, is not guaranteed to match the ordering of the canonical correlations, leading similarly to an inadmissible test.

Nuisance variables

Few authors discussed nuisance variables or confounds in canonical correlation analysis, e.g., Roy (1957); Rao (1969); Timm and Carlson (1976); Lee (1978); Sato et al. (2010). Let the be a matrix of nuisance variables, including an intercept. Partial cca consists of considering nuisance for both and . This is distinct from part cca, which consists of considering a nuisance for either or , but not both. Finally, bipartial cca consists of considering a nuisance for , while considering another set of variables , of size , a nuisance for . In all three cases, such nuisance variables can be regressed out from the respective set of variables of interest, then the respective residuals subjected to cca (Table 1). In the parametric case, inference can proceed using the distribution of or (Equations (3), (4))) using for partial or part, and for bipartial cca (Timm and Carlson, 1976; Lee, 1978).

Table 1

Taxonomy of canonical correlation analysis with respect to nuisance variables. In all cases, the aim is to find linear combinations of variables in the left and in the right sets, such that each combination from one set is maximally correlated with a corresponding combination from the other, but uncorrelated with all other combinations from either set.

Name	Left set	Right set
cca (‘‘full’’, no nuisance)	Y	X
Partial cca	R_ZY	R_ZX
Part cca	R_ZY	X
Bipartial cca	R_ZY	R_WX

is a residual forming matrix that considers the nuisance variables in , and is computed as , where the symbol + represents a pseudo-inverse. is computed similarly, considering the nuisance variables in . The choice which set is on left or right side is arbitrary.

Permutation inference, however, requires further considerations, otherwise, as shown in Section 4, results will be invalid. Consider first the case without nuisance variables. Let be the horizontal concatenation of the two sets of variables whose association is being investigated. Both and occupy an N-dimensional space and, therefore, so does . A random permutation of the rows of either of the two sets of variables will not affect their dimensionalities. For example continues to occupy the same N-dimensional space as .

However, residualisation changes this scenario. Let be the residual forming matrix associated with the nuisance variables , with the symbol + representing the Moore–Penrose pseudo-inverse. has the following interesting properties: (symmetry) and (idempotency), both of which will be exploited later. In partial cca, can be regressed out from and by computing and . Let be the concatenation of the residualised sets and with respect to . While occupies an N-dimensional space, occupies a smaller one; its dimensions are, at most, of a size given by the rank of , which is assuming and are of full rank. The same holds for and and, therefore, for and .

Permutation affects these relations: while still occupies a space of N dimensions as the unpermuted , , differently than , may now occupy a space with dimensions anywhere between and N, depending on a given random permutation. With fewer effective observations determined by this lower space after residualisation, and the same number of variables, the sample canonical correlations in the unpermuted case are stochastically larger than in the permuted, which in turn leads to an excess of spuriously small p-values. For not occupying the same space as the original, the permuted data are no longer a similar realisation of the unpermuted, thus violating exchangeability, and specifically causing the distribution of the test statistics to be unduly shifted to the left.

Here the following solution is proposed: using the results from Huh and Jhun (2001), let be a semi-orthogonal basis (Abadir and Magnus, 2005, p. 84) for the column space of constructed via, e.g., spectral or Schur decomposition, such that , , where , and . Then cca on and leads to the same solutions as on and . The reason is that, from Section 2.1, , which is the same as , since, as discussed earlier, is symmetric and idempotent, and . In a similar manner, , and likewise, . While pre-multiplication by does not affect the cca results,2 it changes the dependence structure among the rows of the data: occupies an -dimensional space, and so does , for a permutation matrix of size , such that exchangeability holds, thus allowing a valid permutation test.

The treatment of partial cca, as described above, can be seen as a particular case of bipartial cca in which , that is, the set of nuisance variables in both sides is the same. Of course, for bipartial cca proper, this equality does not necessarily hold, and the two sets may be different in different ways: may be entirely orthogonal to , or some or all variables from one set may be fully represented in the other, either directly (e.g., some of the variables present in both sets), or as linear combinations of one set in the other, or it may be that these two sets are simply not orthogonal. The direct strategy of computing and its respective semi-orthogonal matrix leads to difficulties because, unless , the products and will not have the same number of rows: has , whereas has rows, thus preventing the computation of cca.

A more general solution, that accommodates bipartial and, therefore, is a generalisation for all cases of nuisance variables in cca, consists of randomly permuting rows of and/or using, respectively, permutation matrices and of respective sizes and , therefore permuting in the lower dimensional space where and are exchangeable, then, crucially, reestablishing the original number N of rows using the property that the transpose of a semi-orthogonal matrix is the same as its inverse (), to only then perform cca. Therefore, cca is computed using and . Left and right sides will continue to have rank and respectively, will have already been permuted, and will both have N rows. The procedure is fully symmetric in that, when the permutation matrices and are both identity matrices (of sizes and , respectively), which is equivalent to no permutation, the expressions for each side reduce to the residualised data and . The concatenation has the same rank as that of , thus addressing the above problem of the unpermuted test statistic having a different and stochastically dominant distribution over that of the permuted data. Table 2 summarises the proposed solution for all cases, including part cca.

Table 2

Proposed permutation method for the various cases of cca, with respect to nuisance variables.

Name	Left set	Right set
cca (‘‘full’’, no nuisance)	P_YY	X
Partial cca	P_YQ_Z’Y	Q_Z’X
Part cca	Q_ZP_YQ_Z’Y	P_XX
Bipartial cca	Q_ZP_YQ_Z’Y	Q_WP_XQ_W’X

is a semi-orthogonal basis for the column space of , such that , , where , and . is a similarly defined matrix for the column space of , . The bipartial cca case generalizes all others: for ‘‘full’’ cca, , and so, ; for partial cca, ; for part cca, and so, . For full and partial, pre-multiplication by can be omitted since , such that results do not change. Once these simplifications are considered, the general bipartial cca case reduces to the other three as shown in the Table. Full and partial have matching number of rows in both sides, such that only one side needs be permuted; part and bipartial, however, have at the time of the permutation a different number of rows in each side, such that both can be permuted separately through the use of suitably sized permutation matrices and ; is size for full cca, and for the three other cases; is size for full and for part cca, and for the two other cases.

Restricted exchangeability

The above method uses the Huh–Jhun semi-orthogonal matrix applied to cca and leads to a valid permutation test provided that there are no dependencies among the rows of . That is, the method takes into account dependencies introduced by the regression of and/or out from and/or , but not dependencies that might already exist in the data, and which generally preclude direct use of permutation tests. However, structured dependencies, such as those that may exist, for instance, in studies that involve repeated measurements, or for those in which participants do not constitute independent observations, e.g., sib-pairs, as in the Human Connectome Project (hcp; Van Essen et al., 2012), can be treated by allowing only those permutations that respect such dependency structure (Winkler et al., 2015). Unfortunately, the Huh–Jhun semi-orthogonal matrix does not respect such structure, blurring information from observations across blocks, and preventing the definition of a meaningful mapping from the N original observations that define the block structure to the or observations that are ultimately permuted.3

Such mapping, whereby each one of the and rows of, respectively, and corresponds uniquely to one of the N rows of the original data and , can be obtained using a different method, due to Theil (1965, 1968), and reviewed in detail by Magnus and Sinha (2005). Consider first the case of . In the Theil method, that here is adapted for cca, , where the exponent represents the positive definite matrix square root, and is a selection matrix, , that is, an identity matrix from which some rows have been removed. Pre-multiplication of a matrix by a selection matrix deletes specific rows, i.e., the ones that correspond to columns that are all zero in the selection matrix (Fig. 1). The thusly computed are the best linear unbiased residuals with scalar covariance (blus), in that they are unbiased estimates of , where are the (unknown) true errors after the nuisance effects of have been removed from ; contains the variance of interest, which may be shared among linear combinations of variables in both sides; it is an estimate of that is subjected to cca and statistical testing. For partial cca, is the same for both sides; for bipartial cca, similar computations hold for the other side, i.e., . Table 3 summarises the two methods.

Fig. 1

A selection matrix is an identity matrix from which some specific rows have been removed. Pre-multiplication by a selection matrix deletes specific rows (those that correspond to columns that are all zero in the selection matrix).

Table 3

The semi-orthogonal matrix ( and/or , subscripts dropped), discussed in Sections 2.6, 2.7, is not unique. Two principled methods to obtain it are below.

Method	Matrix
Theil (1965)	Q=RS’(SRS’)^−1/2
Huh and Jhun (2001)	QEQ’=R (via svd or Schur)

is the residual-forming matrix ( or , for the respective set of nuisance variables, subscript dropped); since is idempotent, all its eigenvalues (the diagonal elements of ) are equal to 0 or 1. In the Theil method, is a (for ) or (for ) selection matrix; the matrix square root (in the exponent ) is the positive definite solution. In the Huh–Jhun method, after Schur or svd factorisations of are computed, and the R or S columns of that have corresponding zero eigenvalues in the diagonal of are removed, such that computed from the factosisation is reduced from size to or to . At the end of these computations (see Algorithm 2, ‘‘semiortho’’, in the Appendix), for both methods, , , and . Both methods aim at obtaining residuals with a scalar covariance matrix . Theil explictly seeks blus residuals. However, strictly, does not need to be a selection matrix: choose to be (not to be confused with qr decomposition) using computed with the Huh–Jhun approach. Then, following Magnus and Sinha (2005, Theorem 2, p. 42), it can be shown that Huh–Jhun also provides blus residuals.

To construct a permutation procedure for cca that respects the block structure, the Theil method can be used to compute instead of the Huh–Jhun approach. Choose observations to be removed from both sides (for partial cca, since ). Construct the selection matrix of size , define the exchangeability blocks based on observations, compute and using the same for both (for part cca, use the same strategy as for bipartial, replacing for ), residualise (in the blus sense) the input variables by computing and . These have the same number of rows, and the dependencies among these rows is the same for both sides; thus, only one side needs to be subjected to random permutations that respect such existing dependencies. Optionally, after permutation, the number N of observations may be reestablished by pre-multiplication by and . Finally, cca is performed, with observation to the aspects discussed in Sections 2.3, 2.5. A detailed algorithm is presented in Section 2.8.

It remains to be decided how to select the observations to be dropped. In principle, any set could be considered for removal, provided that the removed rows of or form a full rank matrix. Some informed choices, however, could be more powerful than others. One of the conclusions from Winkler et al. (2015) is that the complexity of the dependence structure and the ensuing restrictions on exchangeability leads to reductions in power. Thus, natural candidates for removal are observations that, once removed, cause the overall dependence structure to be simpler. For example, it is sometimes the case that some observations are so uniquely related to all others that there are no other observations like them in the sample. These observations, therefore, cannot be permuted with any other, or perhaps with only a few. Their contribution to hypothesis testing in the permutation case is minimal, and their removal is less likely to affect a decision on rejection of the null hypothesis. Consider for example a design that has many monozygotic, dizygotic, and non-twin pairs of subjects, and that in the sample, there happens to be a single pair of half-siblings. It is well known that, for heritable traits, genetic resemblance depends on the kinship among individuals; half-siblings are expected to have a different degree of statistical dependency among each other compared to each one of the other types of sibships in this sample. Thus, in there being just one such pair, it would be reasonable to prioritise it for exclusion, while keeping others.

General algorithm

A set of steps for permutation inference for cca is described in Algorithm 1. In it, input variables and will have been mean-centered before the algorithm begins, or an intercept will have been included as nuisance variable in both and . is a set containing pairs of permutation matrices indexed by . In this set, the first permutation is always ‘‘no permutation’’, i.e., , such that , for all k. For the cases in which only one side of cca needs be permuted (Table 2), or for the cases in which , or when there are dependencies among the data such that the Theil method is used to construct (Table 3), then can be set as for all j. Details on how is defined in observance to the null hypothesis and respecting structured dependencies among the data have been discussed in Winkler et al. (2014, 2015). In the algorithm, P can be larger, equal, or smaller than Q. Optional input arguments are the matrices with nuisance variables and , and the selection matrix . If is supplied but not , then the algorithm performs part or partial cca, depending on the Boolean argument partial; if both and are supplied, the algorithm performs bipartial cca; if neither is supplied, then ‘‘full’’ cca is performed. If is supplied, then the blus residuals based on Theil are used; otherwise, Huh–Jhun residuals are used. For either of these two cases, the semi-orthogonal matrix is computed using a separate, ancillary function named ‘‘semiortho’’, described in the Appendix.

An initial cca using residualised data is done in line 19; this uses another ancillary function, named ‘‘cca’’, and also described in the Appendix; this function returns three results: the canonical coefficients and , and the canonical correlations . The canonical coefficients are used to compute the canonical variables and , augmented by their orthogonal complement needed to ensure that they span the same space as the variables subjected to this initial cca; the canonical correlations are ignored at this point and not stored (hence the placeholder ‘‘_’’). A counter for each canonical component is initialised as 0.

The core part of the algorithm are the two loops that run over the permutations in and the K canonical components (between lines 25 and 35). At each permutation j, cca is executed K times. In each, the columns of and that precede the current k are removed, such that their respective variances are not allowed to influence the canonical correlations at position k. At each permutation, the K canonical correlations are obtained (the third output from the function ‘‘cca’’) and used to compute the associated test statistic. As shown, Wilks’ statistic, , is used, simplified by the removal of the constant term, which does not affect permutation p-values. For numerical stability, sum of logarithms is favoured over the logarithm of a product (compare line 30 with Equation (3)). For inference using Roy’s statistic, replace the condition for in line 31; this modification alone is sufficient as is permutationally equivalent to . In that case, computations indicated in line 30 are no longer needed and can be removed to save computational time.

Whenever the statistic for the correlation at position k in a given permutation is higher or equal than that for the unpermuted data, the counter is incremented (line 32). After the loop, the counter is converted into a p-value for each k. These simple, uncorrected p-values, however, are not useful. Instead, fwer-adjusted p-values are computed under closure using the cumulative maximum, i.e., the p-value for is the largest (least significant) uncorrected p-value up to position k. The algorithm returns then these adjusted p-values, which can be compared to a predefined test level α to establish significance. Note that α itself is never used in the algorithm.

Permutation inference for CCA.

As presented, the algorithm does not cover dimensionality reduction or any penalty to enforce sparse solutions for cca. Dimensionality reduction using methods such as principal component analysis (pca), if included, would be performed after residualisation, but before cca. Thus, in the algorithm, pca, if executed, would be done between lines 18 and 19. As for the many forms of sparse or penalised cca (Nielsen, 2002; Waaijenborg and Zwinderman, 2007; Wiesel et al., 2008; Parkhomenko et al., 2007, 2009; Witten et al., 2009; Soneson et al., 2010; Hardoon and Shawe-Taylor, 2011; Gao et al., 2017; Ma and Li, 2018; Tan et al., 2018), in principle these can be incorporated into the algorithm through the replacement of the classical cca in lines 19 and 29 for one of these methods.

Evaluation methods

In this section we describe the synthetic data and methods used to investigate error rates and power under the different choices for the various aspects presented in Section 2 at each stage of a permutation test for cca, providing empirical evidence for the approach proposed. An overview of these aspects and choices at each stage is shown in Table 4. For each case, we use a series of simulation scenarios: each consists of a set of synthetic variables constructed using random values drawn from a normal or a non-normal (kurtotic or binary) probability distribution, sometimes with or without dimensionality reduction using principal components analysis (pca; Hotelling, 1933; Jolliffe, 2002), sometimes with or without signal, and sometimes with or without nuisance variables. We also consider cases with large sample sizes and large number of variables. An overview of these scenarios (there are twenty of them) is in Table 5.

Table 4

A valid permutation test for cca needs to consider several computational and statistical aspects, which can be seen as steps in a procedure; each of them may be addressed using different strategies, some of which are studied below. The recommendations for the different options for each case (column “Use”) are based on theoretical grounds shown in the sections as indicated in the table, and verified empirically using synthetic data (Section 4).

Step	Possible strategies studied	Use	Theory	Scenarios
Estimation of the canonical components	(a) All in a single step.	Image 6	Section 2.3.	i–vi.
	(b) Stepwise; variance already explained removed.	Image 2
Inclusion of the complement of the canonical coefficients	(a) Null space not included.	Image 6	Section 2.3.	i–vi.
	(b) Null space included.	Image 2
Correction for multiple testing	(a) Uncorrected, simple p-values, [pk]unc.	Image 6	Section 2.5.	i–vi, xviii.
	(b) Corrected, cumulative maximum, [pk]clo.	Image 2
	(c) Corrected, distribution of the maximum, [pk]max.	Image 5
Treatment of nuisance variables	(a) Simple residualisation (Q = I).	Image 6	Sections 2.6, 2.7.	vii–xviii.
	(b) Huh–Jhun method.	Image 2
	(c) Theil method.	Image 2
Choice of the test statistic	(a) Wilks’ λk.	Image 2	Sections 2.2, 2.4.	xvii, xviii.
	(b) Roy’s largest root, θk.	Image 2

Can or should be used.

Can but should not be used.

Cannot or should not be used.

Table 5

Simulation scenarios.

Scenarios		N	P	Q	R	S	#(pca)	Distribution	Signals	#Perms.	#Reps.
Without nuisance	i	100	16	20	0	0	–	normal	–	2000	2000
	ii	100	16	20	0	0	10	normal	–	2000	2000
	iii	100	16	20	0	0	–	kurtotic	–	2000	2000
	iv	100	16	20	0	0	10	kurtotic	–	2000	2000
	v	100	16	20	0	0	–	binary	–	2000	2000
	vi	100	16	20	0	0	10	binary	–	2000	2000
Partial cca	vii	100	16	20	15	R	–	normal	–	2000	2000
	viii	100	16	20	15	R	10	normal	–	2000	2000
	ix	100	16	20	15	R	–	kurtotic	–	2000	2000
	x	100	16	20	15	R	10	kurtotic	–	2000	2000
	xi	100	16	20	15	R	–	binary	–	2000	2000
	xii	100	16	20	15	R	10	binary	–	2000	2000
Bipartial cca	xiii	100	16	20	15	15	–	normal	–	2000	2000
	xiv	100	16	20	15	15	10	normal	–	2000	2000
Larger samples	xv	∗	16	20	20	R	–	normal	–	1000	1000
	xvi	∗	16	20	20	R	10	normal	–	1000	1000
With signal	xvii	100	16	20	0	0	–	normal	sparse	2000	2000
	xviii	100	16	20	0	0	–	normal	dense	2000	2000

∗ For scenarios xv and xvi, sample size varied, . In the table, R and S refer to the number of nuisance variables other than the intercept, which is always included (so the number of nuisance variables in left and right sides for all the simulation scenarios was always, respectively and ). For partial cca, the number of nuisance variables on one side is always the same as in the other, i.e., , but that does not have to be for bipartial cca, even though here the same size was used. The case with larger samples was used for investigation of partial cca.

We start by investigating aspects related to the estimation of the canonical components at each permutation. Specifically, we consider (a) a one-step estimation of all canonical components, from 1 to K, versus (b) sequential estimation that removes, for the k-th canonical component in a given permutation, the variance already explained by the previous ones, as described in Section 2.3. With respect to the inclusion of the complement of the canonical coefficients, we consider (a) without the inclusion of the null space of the canonical coefficients, versus (b) with its inclusion so as to ensure that all variance from the original data not explained in the initial cca is considered in the estimation at every permutation, as described in Section 2.3. With respect to multiple testing, we consider the following strategies: (a) simple, uncorrected p-values, , (b) corrected under closure, , and (c) corrected using the distribution of the maximum statistic ; both and offer fwer control, as discussed in Section 2.5. Keeping the same notation, we define scenarios i–vi consisting of observations, with variables on the left side () of cca and variables on the right side () (the procedure is symmetric; the choice of sides is arbitrary and does not affect results); for these six scenarios, data are drawn from one of three possible distributions: a normal distribution with zero mean and unit variance, a Student’s t distribution with variable degrees of freedom (kurtotic), or a Bernoulli distribution with parameter (binary). Analyses with and without dimensionality reduction to 10 variables using pca are considered. The number of permutations used to compute p-values was set as , with 2000 realisations (repetitions), thus allowing the computation of error rates.

We then turn our attention to aspects related to nuisance and residualisation discussed in Sections 2.6, 2.7. We consider (a) simple residualisation, (b) residualisation using the Huh–Jhun method, and (c) residualisation using the Theil method. For this purpose, scenarios vii–xvii are constructed similarly as i–vi, except that a third set of variables is used as nuisance for partial cca, whereas two other scenarios, xiii and xiv, a fourth set of variables of variables is used as nuisance for bipartial cca.

The impact of ignoring, in samples substantially larger than the number of variables, the dependencies introduced by the residualisation of both sides of cca is studied with scenarios xv and xvi, which consider samples progressively larger, , while keeping the other parameters similar as in scenarios vii and viii. Finally, we briefly investigate power and the choice of the test statistic: we consider (a) Wilks’ statistic (), as well as (b) Roy’s largest root (), as discussed in Sections 2.2, 2.4. We define scenarios xvii and xviii similarly as i, this time including a strong, true signal in one canonical component, thus named ‘‘sparse’’, or multiple, weaker signals shared across multiple (half of the smaller set, thus, ‘‘dense’’). For all scenarios, an intercept is always included as nuisance variable in both sides such that the actual number of nuisance variables is and for each side, respectively. To report confidence intervals (95%), the Wilson (1927) method is used.

Results

In the results below, Sections 4.1, 4.2 establish empirically that with an estimation method (i) that includes the null space of the canonical coefficients, (ii) that finds the canonical correlations in an iterative manner, and (iii) that after computing p-values through a closed testing procedure, the error rates are controlled. The subsequent results, from Section 4.3 onwards, consider only this valid approach.

Estimation strategies

Not including the complement of the canonical coefficients (null space) caused error rates to be dramatically inflated, well above the expected test level (5%), regardless of whether the estimation used the single step or the stepwise procedure, and regardless of any of the multiple testing correction strategies discussed; these results are shown in Table 6.

Table 6

Observed per comparison error rate (%) and 95% confidence intervals for the first 6 canonical correlations in scenario i, assessed using the Wilks’ statistic and three different multiple testing correction methods; the observed familywise error rate (fwer) for each case is also shown. Valid methods should have a fwer close to the nominal 5%, and pcer close to the nominal ; see Sections 4.1, 4.2 for details.

	Null space not included		Null space included
	Single step	Stepwise	Single step	Stepwise
(a) Uncorrected, simple p-values,[p_k]_unc
k=1	91.35 (90.04–92.50)	91.35 (90.04–92.50)	4.70 (3.86–5.72)	4.70 (3.86–5.72)
k=2	93.50 (92.33–94.50)	60.40 (58.24–62.52)	4.60 (3.77–5.61)	0.25 (0.11–0.58)
k=3	94.70 (93.63–95.60)	26.70 (24.81–28.68)	4.60 (3.77–5.61)	0.00 (0.00–0.19)
k=4	95.55 (94.56–96.37)	7.25 (6.19–8.47)	4.85 (3.99–5.88)	0.00 (0.00–0.19)
k=5	96.10 (95.16–96.86)	1.45 (1.01–2.07)	4.40 (3.59–5.39)	0.00 (0.00–0.19)
k=6	96.75 (95.88–97.44)	0.25 (0.11–0.58)	4.30 (3.50–5.28)	0.00 (0.00–0.19)
fwer	99.90 (99.64–99.97)	91.35 (90.04–92.50)	18.30 (16.67–20.05)	4.70 (3.86–5.72)
(b) Corrected, cumulative maximum,[p_k]_clo
k=1	91.35 (90.04–92.50)	91.35 (90.04–92.50)	4.70 (3.86–5.72)	4.70 (3.86–5.72)
k=2	90.35 (88.98–91.57)	60.40 (58.24–62.52)	3.40 (2.69–4.29)	0.25 (0.11–0.58)
k=3	89.80 (88.40–91.05)	26.70 (24.81–28.68)	2.75 (2.12–3.56)	0.00 (0.00–0.19)
k=4	89.55 (88.13–90.82)	7.25 (6.19–8.47)	2.40 (1.81–3.17)	0.00 (0.00–0.19)
k=5	89.30 (87.87–90.58)	1.45 (1.01–2.07)	1.75 (1.26–2.42)	0.00 (0.00–0.19)
k=6	88.85 (87.40–90.16)	0.25 (0.11–0.58)	1.45 (1.01–2.07)	0.00 (0.00–0.19)
fwer	91.35 (90.04–92.50)	91.35 (90.04–92.50)	4.70 (3.86–5.72)	4.70 (3.86–5.72)
(c) Corrected, distribution of the maximum,[p_k]_max
k=1	91.35 (90.04–92.50)	91.35 (90.04–92.50)	4.70 (3.86–5.72)	4.70 (3.86–5.72)
k=2	7.95 (6.84–9.22)	10.95 (9.66–12.39)	0.00 (0.00–0.19)	0.00 (0.00–0.19)
k=3	0.00 (0.00–0.19)	0.00 (0.00–0.19)	0.00 (0.00–0.19)	0.00 (0.00–0.19)
k=4	0.00 (0.00–0.19)	0.00 (0.00–0.19)	0.00 (0.00–0.19)	0.00 (0.00–0.19)
k=5	0.00 (0.00–0.19)	0.00 (0.00–0.19)	0.00 (0.00–0.19)	0.00 (0.00–0.19)
k=6	0.00 (0.00–0.19)	0.00 (0.00–0.19)	0.00 (0.00–0.19)	0.00 (0.00–0.19)
fwer	91.35 (90.04–92.50)	91.35 (90.04–92.50)	4.70 (3.86–5.72)	4.70 (3.86–5.72)

Using the Roy’s statistic led to similar results as with Wilks (not shown). Dimensionality reduction with pca led to similar results for the case in which the null space is included (not shown). For the case in which the null space is not included, results are not comparable with the ones above because, after pca, in the simulations, such that there is no null space to be considered as the matrices with canonical coefficients in both sides are then square.

Even when the null space of the canonical coefficients was included, a single step procedure was never satisfactory. To understand this, consider the following consequence of the theory presented in Sections 2.2, 2.5: for a valid, exact test in cca, the expected error rate for each , i.e., the per comparison error rate (pcer; Hochberg and Tamhane, 1987) is α for , but for it is , since the null can only be rejected if the previous one has also been declared significant at α. More generally, the pcer for a valid test is for the k-th test, i.e., for the k-th canonical correlation. If the test level is set at 5%, then the pcer is 5% for , 0.25% for , 0.0125% for , and so forth. Error rates above this expectation render the test invalid; below render it conservative. In the simulations, a single step procedure never led to an exact test, with or without consideration to multiple testing, as shown in Table 6.

Multiple testing

As with the pcer, it is worth mentioning what the expected fwer for a valid, exact test is. That expectation is the test level itself, i.e., α. Any higher error rate renders a test invalid; lower error rate renders it conservative, though valid. Table 6 shows the fwer for the three different correction methods considered.

If the null space was not included, since the pcer was not controlled, the fwer could not be controlled either (first two columns of the table). If the null space of the canonical coefficients was included (last two columns), even though the single step estimation controlled the pcer, the fwer was not controlled for the simple, uncorrected p-values (third column, upper panel), which is not surprising. It should be emphasised, however, that these simple p-values have another problem: they are not guaranteed to be monotonically related to the respective canonical correlations, such that it is possible that, using these p-values, the null hypothesis could be rejected for some canonical correlation, but retained for another that happens to be larger than the former. The use of such uncorrected, simple p-values, therefore, constitutes a test that is inadmissible. The problem with lack of monotonicity with uncorrected p-values is less severe if estimation is done in a stepwise manner (fourth column, upper panel), but is nonetheless still present, as shown in Fig. 2, and has potential to lead to an excess fwer, even though that did not occur in these simulations.

Fig. 2

Relationship between canonical correlations (horizontal axes) and associated p-values (vertical axes) for 10 realisations of scenario i, considering two estimation methods (single step and stepwise) and three multiple testing correction methods (uncorrected, corrected using the cumulative maximum, and corrected using the distribution of the maximum statistic). The figure complements Table 6 by showing example realisations that average to the error rates shown in the table for the cases in which the null space is included. For simple, uncorrected p-values, the test is inadmissible; for corrected using the distribution of the maximum statistic, the test is overly conservative; single step does not control the familywise error rate.

For the other two correction methods, when the null space of the canonical coefficients was included in the estimation process, fwer was controlled (third and fourth columns of Table 6, middle and lower panels), but there are particularities. Using the distribution of the maximum (lower panel) led to very conservative pcer, for both single step or stepwise estimation, whereas correction with closure led to invalid pcer for single step estimation (third column, middle panel).

The only configuration that led to exact (neither conservative or invalid) control over pcer and fwer, and a monotonic relationship between canonical correlations and associated p-values, is the one in which a stepwise estimation was performed, with the null space of the canonical coefficients included, and with correction using a closed testing procedure (fourth column, middle panel of Table 6). Moreover, the fwer, when controlled using the cumulative maximum or the distribution of the maximum statistic, is guaranteed to match the pcer for : in the former case, any further rejection of the null is conditional on the first one having been rejected; in the latter, the distribution of the maximum coincides with the distribution of the first as the canonical correlations are ranked from largest to smallest.

For partial cca, simple residualisation, even using the above procedure (stepwise estimation, null space included, correction via closure), resulted in the error rates being dramatically inflated, as shown in Table 7. The Huh–Jhun and the Theil residualisation methods, in contrast, resulted in the error rates being controlled at the nominal level, with no excess false positives. For bipartial cca, the problem did not happen in the simulation settings: simple residualisation of both sides by entirely different sets of variables did not cause the error rates to be inflated; yet, using Huh-Jhun or Theil also produced nominal error rates, suggesting that these could be used in any configuration of nuisance variables, regardless of whether those in one side are not independent from those in the other.

Table 7

Observed per comparison error rate (pcer, %) and 95% confidence intervals for the first 6 canonical correlations in scenarios vii (partial cca) and xiii (bipartial cca), for the three different methods considered for treatment of nuisance variables.

	Simple residuals	Huh–Jhun	Theil
(a) Partialcca
k=1 (fwer)	83.85 (82.17–85.40)	5.10 (4.22–6.15)	4.85 (3.99–5.88)
k=2	44.15 (41.99–46.34)	0.30 (0.14–0.65)	0.35 (0.17–0.72)
k=3	12.75 (11.36–14.28)	0.05 (0.01–0.28)	0.00 (0.00–0.19)
k=4	1.75 (1.26–2.42)	0.00 (0.00–0.19)	0.00 (0.00–0.19)
k=5	0.20 (0.08–0.51)	0.00 (0.00–0.19)	0.00 (0.00–0.19)
k=6	0.00 (0.00–0.19)	0.00 (0.00–0.19)	0.00 (0.00–0.19)
(b) Bipartialcca
k=1 (fwer)	5.55 (4.63–6.64)	5.20 (4.31–6.26)	4.45 (3.63–5.44)
k=2	0.10 (0.03–0.36)	0.30 (0.14–0.65)	0.20 (0.08–0.51)
k=3	0.00 (0.00–0.19)	0.00 (0.00–0.19)	0.00 (0.00–0.19)
k=4	0.00 (0.00–0.19)	0.00 (0.00–0.19)	0.00 (0.00–0.19)
k=5	0.00 (0.00–0.19)	0.00 (0.00–0.19)	0.00 (0.00–0.19)
k=6	0.00 (0.00–0.19)	0.00 (0.00–0.19)	0.00 (0.00–0.19)

Estimation included the null space of the canonical coefficients and a stepwise procedure, assessed using the Wilks’ statistic, and corrected using a closed testing procedure (ctp). The ctp guarantees that the familywise error rate (fwer) matches the pcer for the first canonical correlation (i.e., for ). Using the Roy’s statistic led to similar results as with Wilks’; likewise, dimensionality reduction with pca led to similar results (not shown).

Non-normality

Without nuisance variables and with kurtotic data simulated using a Student’s t distribution with a small number of degrees of freedom, , as well as with binary data simulated using a Bernoulli distribution with parameter , error rates were controlled nominally, as shown in Table 8. In partial cca, however, even using the Huh–Jhun method, highly kurtotic data led to excess error rates. In particular, for the simulated data using a Student’s t distribution with degrees of freedom of only , the observed error rate was 14.7%, for a test level of 5%; using the Theil method led to also inflated error rate in this case, with 10.7% (95% confidence interval: 8.28–13.72, not shown in the table). For , error rates were controlled at the nominal level, for both Huh–Jhun (Table 8) and Theil (not shown).

Table 8

Observed per comparison error rate (pcer, %) and 95% confidence intervals for the first canonical correlation in scenarios i–vi (without nuisance) and vi–xii (partial cca, with Huh–Jhun), considering different distributions for the data.

Distribution		Without nuisance	Partial cca
Normal		4.70 (3.86–5.72)	5.15 (4.26–6.21)
Student	ν=2	3.95 (3.18–4.90)	14.70 (13.22–16.32)
	ν=4	5.45 (4.54–6.53)	5.40 (4.49–6.48)
	ν=6	4.15 (3.36–5.12)	5.40 (4.49–6.48)
	ν=8	4.70 (3.86–5.72)	5.00 (4.13–6.04)
	ν=10	3.85 (3.09–4.79)	5.10 (4.22–6.15)
Bernoulli	q=0.2	5.30 (4.40–6.37)	5.30 (4.40–6.37)

: Degrees of freedom of the Student’s t distribution used to simulate data; q: Parameter of the Bernoulli distribution used to simulate data. Estimation used the null space of the canonical coefficients and a stepwise procedure, assessed using the Wilks’ statistic, and corrected using a closed testing procedure (ctp). The ctp guarantees that the familywise error rate (fwer) matches the pcer for the first canonical correlation (i.e., for ). Using the Roy’s statistic led to similar results as with Wilks’; using Theil led to similar results as Huh–Jhun; likewise, dimensionality reduction with pca led to similar results (not shown).

Large samples

Increasing the sample size while keeping the number of variables fixed progressively reduced the amount of errors for the simple residualisation method to treat nuisance variables, as shown in Table 9; the trend was similar with or without dimensionality reduction using pca. The reduction in the error rate as the sample size increased did not affect Huh–Jhun or Theil methods, for which error rates were already controlled even with a relatively smaller sample compared to the number of variables.

Table 9

Observed familywise error rate (fwer, %, that matches the pcer for ) and 95% confidence intervals for scenarios xv and xvi, used to investigate effect of different sample sizes, for the three different methods for dealing with nuisance variables.

N	Without pca			With pca
	Simple residuals	Huh–Jhun	Theil	Simple residuals	Huh–Jhun	Theil
100	96.60 (95.29–97.56)	4.80 (3.64–6.31)	5.00 (3.81–6.53)	59.20 (56.12–62.21)	5.00 (3.81–6.53)	5.10 (3.90–6.64)
200	42.20 (39.17–45.29)	5.50 (4.25–7.09)	5.40 (4.16–6.98)	22.30 (19.83–24.98)	5.00 (3.81–6.53)	4.50 (3.38–5.97)
300	25.10 (22.51–27.88)	5.00 (3.81–6.53)	5.40 (4.16–6.98)	12.20 (10.31–14.37)	5.40 (4.16–6.98)	5.10 (3.90–6.64)
400	18.00 (15.74–20.50)	4.30 (3.21–5.74)	5.10 (3.90–6.64)	11.80 (9.95–13.95)	4.50 (3.38–5.97)	5.70 (4.43–7.31)
500	11.50 (9.67–13.63)	6.10 (4.78–7.76)	4.10 (3.04–5.51)	9.50 (7.83–11.48)	6.80 (5.40–8.53)	5.10 (3.90–6.64)
600	10.80 (9.02–12.88)	5.20 (3.99–6.76)	5.00 (3.81–6.53)	9.30 (7.65–11.26)	4.70 (3.55–6.19)	4.70 (3.55–6.19)
700	11.20 (9.39–13.31)	4.20 (3.12–5.63)	4.40 (3.29–5.86)	7.20 (5.76–8.97)	5.70 (4.43–7.31)	5.50 (4.25–7.09)
800	10.10 (8.38–12.12)	5.50 (4.25–7.09)	4.20 (3.12–5.63)	7.00 (5.58–8.75)	4.80 (3.64–6.31)	5.90 (4.60–7.54)
900	8.40 (6.84–10.28)	5.00 (3.81–6.53)	4.00 (2.95–5.40)	7.70 (6.20–9.52)	5.20 (3.99–6.76)	5.20 (3.99–6.76)
1000	7.70 (6.20–9.52)	4.30 (3.21–5.74)	5.90 (4.60–7.54)	5.00 (3.81–6.53)	6.20 (4.87–7.87)	4.70 (3.55–6.19)

Without pca: , . With pca: . Estimation included the null space of the canonical coefficients and a stepwise procedure, assessed using the Wilks’ statistic, and corrected using a closed testing procedure (ctp). The ctp guarantees that the familywise error rate (fwer) matches the pcer for the first canonical correlation (i.e., for ). Using the Roy’s statistic led to similar results as with Wilks’. The confidence intervals are wider than for other tables because the number of realisations (and also of permutations) was smaller (Table 5).

Dimensionality reduction

Dimensionality reduction with pca did not affect error rates (pcer and fwer) with respect to single step vs. stepwise estimation of canonical coefficients, nor correction for multiple testing, nor method for addressing nuisance variables. That is, these results (not shown) were indistinguishable from those obtained without pca (shown above). Moreover, as the simulations used the same number of principal components for both sides of cca, including or not the null space could not have affected results, as after dimensionality reduction. Using pca did yield higher power to detect effects, for both Wilks’ and Roy’s test statistics (Table 10, next item). This apparent extra power can be attributed to the smaller number of variables after pca, as the principal components that were retained contained most of the simulated signal, which, given the reduced dimensionality of the set of data, could then be detected with higher likelihood.

Table 10

Observed power (%) and 95% confidence intervals for the first canonical correlation in scenarios xvii and xviii, that included a synthetic signal added to either one or half (‘‘sparse’’ or ‘‘dense’’, respectively) of the initial variables, thus captured by only one (sparse case) or multiple (dense case) canonical correlations.

Signals		Without pca		With pca
		Wilks (λ)	Roy (θ)	Wilks (λ)	Roy (θ)
Sparse	k=1	42.10 (39.95–44.28)	57.90 (55.72–60.05)	81.55 (79.79–83.19)	94.75 (93.68–95.64)
Dense	k=1	83.05 (81.34–84.63)	37.05 (34.96–39.19)	95.80 (94.83–96.59)	70.70 (68.67–72.65)
	k=2	42.30 (40.15–44.48)	4.75 (3.90–5.77)	72.05 (70.04–73.97)	24.75 (22.91–26.69)
	k=3	12.75 (11.36–14.28)	0.25 (0.11–0.58)	31.95 (29.94–34.03)	4.30 (3.50–5.28)
	k=4	1.95 (1.43–2.65)	0.00 (0.00–0.19)	6.55 (5.55–7.72)	0.50 (0.27–0.92)
	k=5	0.10 (0.03–0.36)	0.00 (0.00–0.19)	1.00 (0.65–1.54)	0.00 (0.00–0.19)
	k=6	0.05 (0.01–0.28)	0.00 (0.00–0.19)	0.05 (0.01–0.28)	0.00 (0.00–0.19)

Estimation used the null space of the canonical coefficients and a stepwise procedure, assessed using the Wilks’ statistic, and corrected using a closed testing procedure (ctp). The ctp guarantees that the familywise error rate (fwer) matches the pcer for the first canonical correlation (i.e., for ).

The results above, that consider solely the error rates, and are based on results with the Wilks’ statistic (), are essentially the same for Roy’s largest root (; results not shown). That is, results regarding the estimation strategies, multiple testing, nuisance variables, non-normality, behaviour with large samples, and dimensionality reduction with pca, are virtually the same for Wilks’ and Roy’s statistics. In the presence of synthetic signal, however, the two test statistics diverged. Table 10 shows that, with signal spread across multiple canonical components (i.e., ‘‘dense’’), Wilks’ is substantially more powerful than Roy’s statistic. With signal concentrated in just one (the first) canonical variable (i.e., ‘‘sparse’’), the trend reverses, and Roy’s become more powerful than Wilks’.

Discussion

Permutation tests

Compared to univariate, multivariate tests pose the problem of establishing the distributional form for more complicated test statistics; in the parametric case, inference is marred by a set of difficulties: the assumption that all observations are independent and identically distributed following normal theory, the extremely complicated formulas for the density of the canonical correlations, which further depend on the (unknown) population canonical correlations, the sensitivity of asymptotic approximations to departures from assumptions, bias in estimations of parameters, and the validity of these approximations only for particular cases.

Permutation tests address these difficulties in different ways, and their advantages are well known (Ludbrook and Dudley, 1998; Nichols and Holmes, 2002; Good, 2005; Pesarin and Salmaso, 2012): no underlying distributions need be assumed, non-independence and even heteroscedastic variances can be accommodated, non-random samples can be used, and a wide variety of test statistics are allowed. Moreover, all information needed to build the null distribution lie within the data, as opposed to in some idealised population.

These many benefits extend to inference for multivariate methods. In the case of cca, one benefit is immediately obvious: the complicated formulas and charts for the distribution of the canonical correlations can be bypassed completely, thus with no need to appeal to distributional assumptions. In effect, as shown in Section 4.4, even with all variables not following a normal distribution, error rates were still controlled at the nominal level. It should be noted, however, that extremely kurtotic data, such as that generated with a Student’s t distribution with extremely low degrees of freedom, caused results to be invalid in the presence of nuisance variables, even with the Huh–Jhun or Theil methods. Such data, however, are rare (recall that with 2 degrees of freedom, the Student’s t distribution has infinite variance); most applications of cca investigate datasets that have variables with data that have diverse distributional properties.

Yet, although in the univariate case, algorithms for permutation inference tend to be relatively straightforward to implement and do lead to valid results, for cca, the theory presented in the previous sections and the results with synthetic data show that a simple permutation algorithm that does not consider aspects such as a stepwise estimation of the canonical correlations, nor the inclusion of the null space of canonical coefficients when the two sets of variables do not have the same size, or that does not accommodate specific treatment for nuisance variables, or addresses multiplicity respecting the ordering of the canonical correlations, leads to invalid results.

Estimation and multiple testing correction

Results from Sections 4.1, 4.2 show that the estimation method that leads to exact, valid results (neither conservative or invalid) is the one that estimates one canonical correlation one at a time, in a stepwise, iterative manner, that includes the null space of the canonical coefficients when the sets of variables have different sizes (i.e., when ), and that computes adjusted p-values using a closed testing procedure. All alternative approaches led to either invalid or conservative results when considering pcer or the fwer.

It should be emphasised, however, that there are cases in which the naïve permutation method, described at the beginning of Section 2.3 remains valid. The method is valid whenever only the first () canonical component is of interest, and there are no nuisance variables or, if there are nuisance variables, those in the left and right side ( and ) are completely orthogonal (thus, excluding partial cca). Even though the naïve method was not explicitly tested, it is equivalent to the single step method with the null space included, which in the simulations led to an error rate of 4.70% at test level 0.05 (Table 6). The reason why it remains valid is that, if interest is only in the first canonical component, there is no need to perform an initial cca to allow stepwise removal of previous (before the current k components). Moreover, there is no multiple testing to be considered.

The last column of Table 6 may suggest that uncorrected p-values (upper panel) and a ctp (middle) are equivalent for stepwise estimation. They are not, and their differences are manifest in two ways, both previously discussed: first, uncorrected p-values are not monotonically related to the canonical correlations (Fig. 2), and second, fwer has potential to be higher than the pcer for , even though that did not happen in the simulations.

Inference in the presence of nuisance variables

It is sometimes the case that known, spurious variability needs to be taken into account. For example, variables such as age and sex are often considered confounds. Merely regressing out such nuisance variables from all other variables that are subjected to cca, then proceeding to a simple permutation test, leads to inflated error rates and an invalid test, as expected from Section 2.6, and evidenced by the results in Section 4.3. The dependencies among observations introduced through the residualisation renders the data no longer exchangeable.

This inflated error rate, even after multiple testing correction, is the probably the most striking finding of the current study, as the results can be dramatically affected, particularly if the number of nuisance variables is relatively large compared to the sample size, as shown in Section 4.5. Transformations that make residuals exchangeable, through the use of a lower dimensional basis where exchangeability holds, namely, the Huh–Jhun and Theil methods, mitigate the problem, as evidenced by the theory and through the simulations.

Even though both methods led to similarly controlled error rates, they are not equivalent: Huh–Jhun always leads to same canonical components as they would have been obtained from the residualised data, whereas the Theil method can allow for multiple, different solutions depending on the choice of the selection matrix . Theil (1965) suggested that the choice of the observations to be dropped should consider power; here we suggest that the choice of can be based on restrictions on exchangeability: if all original data are freely exchangeable, the Huh–Jhun method is a preferable choice in that it does not require an additional argument that affect the results; however, it does require a Schur or singular value decomposition of the residual-forming matrix, which is a rank-deficient matrix, such that numerical stability should also be a factor for consideration.

For bipartial cca, while error rates were controlled even in the simple residualisation case, it should be noted that and were generated independently in the simulations, such that they were expected to be orthogonal. With real data, possible overlap among columns or linear combinations of columns between and create a case that would lie between the two extremes of partial and bipartial cca. In such case, and given the results for partial cca, error rates are not expected to be controlled with simple residualisation. Huh–Jhun and Theil, being able to deal with the most extreme case of dependencies between and (that is, when the two are the same, which defines partial cca), constitute a general solution to all cases.

Relationship with the GLM

The dangers of residualising both dependent and independent variables in the general linear model (glm) with respect to nuisance variables, then proceeding to a permutation test, as proposed originally by Kennedy (1995) are well known (Anderson and Robinson, 2001). It is not a complete surprise, therefore, that permutation inference for cca would lead to invalid results in similar settings. The original Huh and Jhun (2001) method (see also Kherad-Pajouh and Renaud, 2010) was proposed for the glm as a way to address shortcomings of the Kennedy method in accommodating nuisance variables. Both Kennedy and Huh–Jhun were evaluated by Winkler et al. (2014): among the methods that can be considered for permutation inference in the glm, Huh–Jhun is the only that cannot be used directly with exchangeability blocks, as the reduction to a lower dimensional space does not respect the block structure. The solution proposed here for permutation inference for cca in the presence of exchangeability blocks, which uses the Theil method, is expected to solve the same problem also for the glm, i.e., as a replacement for Huh–Jhun in cases where the data have a block dependence structure, as it does for freely exchangeable data (Ridgway, 2009).

As in the univariate case, permutation tests in the presence of nuisance variables are approximate. Their exactness is in the sense that, under the null hypothesis, the probability of finding a p-value smaller or equal to the test level is the test level itself. Such tests are not perfectly exact as the true relationship between the nuisance variables and the variables of interest are not known and needs to be estimated. Even in the absence of nuisance variables, however, permutation tests that use only a fraction of the total number of possible permutations are also approximate, for not covering the whole permutation space (the number of potential permutations tends to be very large, and grows very rapidly with increases in sample size). The same holds for other resampling methods that do not use all possible rearrangements of the data. Regardless of the reason why the tests are approximate, results are known to converge asymptotically to the true p-values.

Among the two test statistics considered, Wilks’ () tends to be more powerful than Roy’s () for effects that span multiple canonical components; the converse holds for signals concentrated in only a few of the canonical components, i.e., when many of the canonical variables are zero; in these cases, Roy’s tend to be more powerful than Wilks’, as shown in Section 4.7. The respective formulas (Equations (3), (4))) give insight on why that is the case: Roy’s statistic is invariant to canonical correlations other than the first (largest), whereas Wilks’ pool information across all correlations; past simulations, reviewed by Johnstone and Nadler (2017), corroborate to the finding.

The use of these two statistics for any canonical correlation other than the first (i.e., for ) is possible in the proposed iterative procedure because, for the current position k being tested, all the variance associated with the previous canonical components at positions will have already been removed from the model (Sections 2.3, 2.8), such that the largest canonical correlation (Roy’s statistic) is the current one being tested; for Wilks’, the procedure holds because these earlier canonical correlations are not marked as zero; instead, they are ignored altogether when the statistic is computed, as if the previous canonical components have never existed.

Wilks’ lambda and Roy’s largest root are not the only possible statistics that can be considered for cca, and permutation tests allow the use of yet others. Some, such as Hotelling–Lawley and Pillai–Bartlett, were considered by Friederichs and Hense (2003). Using simulations and Monte Carlo results, the authors found that parametric distributions of these classical multivariate statistics were accurate, and could be obtained quickly at low computational cost; it should be noted, however, that the study used normally distributed simulated data, in which case parametric assumptions are known to hold.

Relationship with previous studies

While a number of studies have used permutation tests with cca, not many investigated the performance of these tests. Nandy and Cordes (2003) proposed a non-parametric strategy for inference with cca for the investigation of task-based fmri time series: the method uses a resting-state (no task-related activity) dataset to build the null distribution; as resampling time series can be challenging due to temporal autocorrelation, the null distribution uses multiple voxels selected far apart from each other so as to also avoid issues with spatial correlation. The approach differs from the one presented here in that it uses subject-level time series (as opposed to between-subject analyses), is specific to brain imaging (the proposed method is general) and does a resampling method that shares similarities with, yet is not the same as permutation. Eklund et al. (2011) specifically used permutation tests for cca with fmri time series whitened with a combination of methods to allow permutation; the authors demonstrated that permutation tests for both the glm and cca could be greatly accelerated through the use of graphics processing units (gpus).

Kazi-Aoual et al. (1995) proposed analytical formulas for the first three moments of the permutation distribution of Pillai’s trace for cca; these moments can be used to fit a Pearson (1895) type iii distribution, from which p-values can be obtained. Legendre et al. (2011) studied parametric and permutation tests for redundancy analysis (rda; Rao, 1964) and for canonical correspondence analysis (referred to also by the acronym ‘‘cca’’; ter Braak, 1986); the authors found that a simultaneous test of all canonical eigenvalues for the respective axes (eigenvectors of predicted response variables in a linear model) in rda, despite simple, is not valid, whereas a marginal test on each eigenvalue, as well as a ‘‘forward’’ test in which previously tested canonical axes are added to a matrix of nuisance variables, performs well, even if conservatively for axes other than the first. Yoo et al. (2012) investigated the relationship between cca and regression, proposing the use of permutations and studying cases without nuisance variables; in the method, for the k-th canonical correlation, variance not already explained by canonical variables in one of the sides is permuted, whereas the other variables remain fixed. Turgeon et al. (2018) considered using a small number of permutations for cca, recording of the empirical distribution function, then using it to estimate the parameters of a Tracy–Widom distribution (Tracy and Widom, 1996; Johnstone, 2008) for cases in which the number of observations is smaller than the number of variables in either or ; the distribution is then used to obtain p-values; data are assumed to follow a normal distribution, and inference is for the largest canonical correlation.

Permutation tests for the method of partial least squares (pls; Tucker, 1958; Wold et al., 1983; McIntosh et al., 1996; McIntosh and Lobaugh, 2004) have been considered. For example, Chin and Dibbern (2010) and Sarstedt et al. (2011) used a permutation test to investigate how differences in the strength of association between variables (magnitude of estimates) further differed between two or more groups. These would be equivalent to, in the context of cca, testing whether canonical correlations obtained across different groups would differ. Le Floch et al. (2012) investigated strategies for dimensionality reduction and regularisation for imaging and genetic data, whereas Grellmann et al. (2015) compared the performance of variants of cca and pls for similar problems. Both studies used direct permutation of the data, and were mostly focused on the relative performance of the different methods, offering no specific treatment of nuisance variables or the other aspects considered here, and which concern validity. Monteiro et al. (2016) investigated a strategy for sparse pls and sparse cca in which data are split into training and hold out, and inference uses permutation of the training data, with coefficients applied to test data, in which measurements of association are computed.

The current paper therefore fills a substantial knowledge gap, whereby not many studies considered at all the validity of permutation inference with cca, but those that did approach the topic were not sufficiently general; none covered the topics discussed here. Moreover, in principle, the method as proposed can be used with subject-level fmri or other timeseries data provided that whitening has been successful in removing temporal dependencies. Additionally, given the conceptual similarity between pls and cca, it is possible that permutation inference for pls would require similar strategies as described in Section 2, particularly in the presence of nuisance variables and re-use of variance already explained. Whether that is the case, it is a question that remains open for future investigation.

Recommendations

Given the above results, the main recommendations for permutation inference for cca can be summarised as follows:

When studying a given k-th canonical variable or canonical correlation, , remove the effects of the previous ones, i.e., the variance from one set that has already been explained by the other, as represented by the earlier canonical variables. These effects are surely significant (regardless of the test level), otherwise the current canonical variable or correlation would not be under consideration. Ignoring the earlier ones causes the error rates be inflated (empirical evidence provided in Section 4.1).

For sets of variables with different sizes (i.e., ), ensure that the variability not represented by the canonical variables produced at the first permutation is considered in all and every permutation. That is, include the null space of the canonical coefficients when computing the variables subjected to permutation. Not including the null space leads to excess false positives (empirical evidence provided in Section 4.1).

Do not use simple p-values for inference, and make sure that a closed testing procedure is used. Using simple, uncorrected p-values has two negative consequences: (i) both the pcer and fwer are inflated, and (ii) since simple p-values are not guaranteed to be monotonically related to the canonical correlations, the resulting test is inadmissible (empirical evidence provided in Section 4.2).

For the same reason, do not use fdr to correct for multiple testing after using simple p-values: while the p-values themselves satisfy the requirements of fdr, they lead to an inadmissible test even after correction, leading to non-sensical results whereby a stronger canonical correlation may be less significant than a weaker one (empirical evidence provided in Section 4.2).

While valid, inference using the distribution of the maximum statistic across canonical correlations leads do conservative results, except for the first canonical correlation (empirical evidence provided in Section 4.2).

If regressing out nuisance variables from both sets of variables subjected to cca, make sure that the residuals are transformed to be exchangeable, e.g., with the Huh–Jhun or Theil methods, then permuted accordingly. Failure to observe this recommendation leads to excess false positives, particularly when the number of nuisance variables is a large fraction of the sample size (empirical evidence provided in Sections 4.3, 4.5).

All these recommendations are integrated into Algorithm 1.

Conclusion

As evidenced by the theory and simulations in the previous sections, a simple permutation procedure leads to invalid results: (i) simple p-values are not admissible for inference in cca, lead to excess pcer and fwer, and cannot be corrected using generic methods based on p-values such as fdr; (ii) ignoring the variability already explained by previous canonical variables leads to inflated error rates for all canonical correlations except for the first; (iii) regression of the same set of nuisance variables from both sides of cca without further consideration leads to inflated error rates; and (iv) the classical method for multiple testing correction, that uses the distribution of the maximum statistic, leads to conservative results. The use of a stepwise estimation procedure, transformation of the residuals to a lower dimensional basis where exchangeability holds, and correction for multiple testing via closure, ensures the validity of permutation inference for cca.