Correlated Components Analysis - Extracting Reliable Dimensions in Multivariate Data
Lucas C. Parra Stefan Haufe Jacek P. Dmochowski
stat.ML cs.LG


[
    [
    December 30, 2023
=====================




How does one find dimensions in multivariate data that are reliably expressed across repetitions? For example, in a brain imaging study one may want to identify combinations of neural signals that are reliably expressed across multiple trials or subjects.  For a behavioral assessment with multiple ratings, one may want to identify an aggregate score that is reliably reproduced across raters.  Correlated Components Analysis (CorrCA) addresses this problem by identifying components that are maximally correlated between repetitions (e.g. trials, subjects, raters). Here we formalize this as the maximization of the ratio of between-repetition to within-repetition covariance. We show that this criterion maximizes repeat-reliability, defined as mean over variance across repeats, and that it leads to CorrCA or to multi-set Canonical Correlation Analysis, depending on the constraints. Surprisingly, we also find that CorrCA is equivalent to Linear Discriminant Analysis for zero-mean signals, which provides an unexpected link between classic concepts of multivariate analysis. We present an exact parametric test of statistical significance based on the F-statistic for normally distributed independent samples, and present and validate shuffle statistics for the case of dependent samples. Regularization and extension to non-linear mappings using kernels are also presented. The algorithms are demonstrated on a series of data analysis applications, and we provide all code and data required to reproduce the results.







§ INTRODUCTION


Consider the following scenario: A group of people are shown a movie while their neural activity is recorded.  Guided by the assumption that the movie evokes similar brain responses in different subjects, the goal is to identify movie-related brain activity that is common across individuals  <cit.>. Because activity is distributed across multiple sensors, we would like to find linear combinations of sensors that have a similar time course in all subjects. In other words, we would like to identify underlying components or factors that have high inter-subject correlation.

In a different setting, the task may be to assess motor skills. A clinician observes a group of individuals performing various tasks and provides a rating for each person and task  <cit.>. Typically, ratings are combined across tasks so that each person obtains a single aggregate performance score. But how should individual task ratings be combined? Traditionally, this is done by simple averaging of related scores. We propose instead a data-driven approach that aims to make the aggregate scores as consistent as possible across different raters. This means that we want to combine scores so as to maximize inter-rater reliability. 

What is common in these two scenarios is the goal of identifying directions in multivariate data sets that are reliably reproduced across repetitions.  Correlated Components Analysis (CorrCA) accomplishes this goal. CorrCA was originally developed in the context of neuroimaging studies to extract similar activity in multiple subjects <cit.>, and was re-developed independently to extract reliable brain responses across repeated trials in a single subject <cit.>. Here we will show that the technique can also be used to identify aggregate ratings with high inter-rater reliability. 

More generally, CorrCA is applicable whenever a data volume is available with dimensions T × D × N, where N are repeated measures with the same sensors. The method identifies directions in the D-dimensional space along which the data maximally correlate between N repeated measurements, with correlation measured across T samples. This is similar to Principal Components Analysis <cit.> in that CorrCA finds a set of D-dimensional projection vectors that linearly decompose a data set (of size T × D).  Instead of capturing dimensions with maximum variance within a single data set as in PCA, CorrCA captures the dimensions with maximum correlation between the N data sets. As we will show here, the optimality criterion of CorrCA can also be used to derive Canonical Correlation Analysis <cit.>, including multi-set CCA <cit.>. The methods differs in that CorrCA uses a single set of linear projections for all data sets assuming they repeated measures are all taken in the same space, whereas CCA and MCCA yield a different set of projection vectors for each data-set. MCCA is well established for the purpose of identifying similar brain activity across subjects <cit.>. Here we will show for brain signals and rating data that CorrCA can achieve better performance than MCCA by virtue of the smaller number of free parameters. 


The main contribution of this paper is to provide a didactic yet formal characterization of CorrCA and with this, establish a few novel theoretical results. In particular, we show that maximizing the correlation between repeated measurements is indeed equivalent to maximizing the reliability of the projected data. By “reliability” we mean the average over repetitions, divided by the variance over repetitions, which is also a sensible definition for signal-to-noise ratio. We establish a direct link to MCCA, as well as Linear Discriminant Analysis (LDA). CorrCA and MCCA both maximize inter-subject correlation, defined here the ratio of between- to within-subjects covariance, while LDA  maximizes the between-class over the within-class variance. Surprisingly, LDA and CorrCA yield the same result for zero-mean signals. This direct equivalence allows us to formulate a strict statistical test for significance based on Fisher's variance ratio (F-statistic). The equivalence to LDA suggests a straightforward extension of CorrCA to non-linear mappings using the well-known “kernel” approach. We also provide a validation of CorrCA for identifying shared dimensions using simulated data and validate the proposed tests of statistical significance. Another important contribution of this paper is to formalize inter-subject correlation (ISC) as a correlation metric which applies to more than two signals. ISC differs from conventional definitions of correlation, yet seamlessly integrates into the frameworks of Linear Discriminant Analysis, Canonical Correlation Analysis, and the F-statistic.

Before we present the mathematical treatment, we begin with a simple two-dimensional example to illustrate the basic concept of identifying reliable signal components that are contained in a shared linear dimension. 






 §.§ Extracting shared dimensions — an illustrative example


For this illustrative example we consider a simulation of electrical brain signals recorded from multiple subjects. Assume that we have two sensors, x_1(t) and x_2(t), which record temporal signals from each of two subjects, as shown in Figure <ref> A & B. The figure caption details the generation of the data; briefly, there are two “sources” of signal: one has a common waveform for both subjects, while a second source has a subject-dependent waveform. The common waveform indroduces a correlation between subjects at both sensors, although the exact shared waveform is not immediately apparent in the individual time courses shown in panels A & B). The source directions become clear from a scatter plot of the same data (panel C, blue and red arrows). Here, for each time point t we have connected the samples of both subjects. Evidently, at any point in time, the signals of the two subjects lie along a line (in the direction of the source not shared by the subjects, blue arrow), but the lines are shifted in parallel (in direction of the common source, red arrow). 
The locations along the gray lines are determined by the subject-dependent waveform, while the shift of the parallel line from one time instance to the next is determined by the common waveform. 
It is important to note that while both directions are shared by the two subjects, only one of these dimensions expresses a common waveform. 
Now we select a linear transformation, V, to map this data onto a component space, y(t)= V^⊤ x(t),
such that the component signals are maximally correlated between subjects. With this choice, the time course for both subjects in the first component dimension, y_1(t), are now perfectly correlated (panel D). Thus, we have extracted the common source of variation in these two subjects, whereas the second component dimension, y_2(t), now captures a waveform that is not necessarily shared between subjects (panel E). As a result, the gray lines in the scatter plot are now vertical (panel F), separating the waveforms into reproducible signal (y_1) and non-reproducible noise dimensions (y_2). Evidently, V captured the source directions in these data (blue and red arrows). The objective of Correlated Components Analysis is to identify the unknown projection matrix V by analyzing time courses x(t). The fundamental assumption of CorrCA is that there is a reproducible signal and a non-reproducible noise, and that the directions of the reproducible signal are shared between subjects. 

Readers wishing to omit the theoretical treatment below may proceed directly to Section <ref>, which presents several applications of CorrCA to real-world data sets.  



§ THEORY




 §.§ Problem statement: Maximal inter-subject correlation



The objective of CorrCA is to find dimensions in multivariate data that maximize correlation between repetitions. Denote the D-dimensional data by x_i^l ∈ℝ^D, where i=1, … ,T enumerates exemplars, and l=1, …, N represents repeated renditions of these data. In the case of brain signals, the exemplars represent time points, so that i indexes a sequence in time of length T, and l indexes N repeated recordings (subjects or trials). In the case of a behavioral assessment instrument, i may index individuals receiving a rating, while l indexes multiple raters, or repeated ratings given to the same individuals by a single rater. In both instances, the observations are D-dimensional (D sensors that record brain signals, or D items rated in the assessment instrument). For simplicity, in the following we will employ the terms “time” and “subjects” in reference to i, such that correlation is measured between subjects and across time. When appropriate, we will rephrase the resulting expressions for the case of assessing inter-rater correlation or test-retest correlation. 

The goal is now to identify a linear combination of these D measures, defined by a projection vector v∈ℝ^D:


    y^l_i =  v^⊤ x^l_i ,
 

such that the correlation between N subjects (repetitions) is maximized.  We define inter-subject correlation (ISC) as the ratio of the between-subject covariance, r_B, and the within-subject covariance, r_W: 


    ρ = 1/N-1r_B/r_W .


The between-subject covariance is a sum over all pairs of subjects, and within-subject covariance is a sum over all subjects: 


    r_B    =   ∑_l=1^N ∑_k=1, k ≠ l^N ∑_i=1^T   (y_i^l - y̅_*^l)(y_i^k - y̅_*^k) 
     , 
    
    r_W    =   ∑_l=1^N ∑_i=1^T (y_i^l - y̅_*^l) (y_i^l - y̅_*^l)
     ,


where, y̅^l_* =  1/T∑_i=1^T y_i^l, is the sample-mean of the projected data for subject l. The factor (N-1)^-1 in definition (<ref>) is required so that correlation is normalized to ρ≤ 1 (see Appendix <ref>). To simplify the notation, in our “variance” and “covariance” definitions we generally omit common scaling factors, in this case we omitted (T-1)^-1N^-1.

While we refer to ρ as the inter-subject correlation (ISC), it could also be termed the inter-rater correlation or inter-repeat correlation (IRC) depending on the context. In Appendix <ref> we discuss this definition of correlation and its relation to other measures of correlation, such as intra-class correlation or Pearson's correlation coefficient. Note that r_W is also the variance around the subject-mean.  


By inserting Eq. (<ref>) into definition (<ref>) it follows readily that


    ρ = 1/N-1 v^⊤ R_B  v/ v^⊤ R_W  v ,


where R_B and R_W are the between-subject and within-subject covariance matrices of x_i^k defined analogously to (<ref>)-(<ref>):


    R_B    =   ∑_i=1^T ∑_l=1^N ∑_k=1, k ≠ l^N 
    ( x_i^l - x̅_*^l)( x_i^k - x̅_*^k) ^⊤ ,
    
     R_W    =   ∑_i=1^T ∑_l=1^N 
    ( x_i^l - x̅_*^l)( x_i^l - x̅_*^l)^⊤ .


Here x̅_*^l = 1/T∑_i=1^T  x_i^l is the sample-mean vector for subject l. The projection vector v that maximizes ρ can be found as the solution of  ∂ρ / ∂ v^⊤ = 0, which yields: 


    R_B  v = ρ  R_W  v   .


Assuming that R_W is invertible, v is an eigenvector of R_W^-1 R_B with eigenvalue ρ. Maximal ISC is achieved by projecting the data onto the eigenvector of R_W^-1 R_B with largest eigenvalue. We refer to this principal eigenvector as v_1. In addition to this projection with maximal ISC, there may be additional dimensions that capture significant ISC in the data. Our goal is to find all such projection vectors, v_d, d=1 … D. Thus, a natural objective is to maximize ISC summed over all components. However, we want each component to capture a different, independent aspect of the data. We therefore require components to be uncorrelated, similar to the constraints enforced by Independent Components Analysis or Canonical Correlation Analysis. (The alternative of requiring orthogonal projection vectors, as in PCA is discussed in Appendix <ref>.)  Thus, in total the goal is to maximize the following cost function with respect to V = [ v_1, … v_K]:


    J( V) = ∑_d=1^Dρ_d = ∑_d=1^D v_d^⊤ R_B  v_d/ v_d^⊤ R_W  v_d


subject to the constraints that components are mutually uncorrelated, i.e., v_d^⊤ R_W  v_c =0, for all c  d. The solution to this problem is given by the  eigenvector matrix V that satisfies the eigenvalue equation (see Appendix <ref>):


    R_B  V =   R_W  VΛ,


where Λ is a diagonal matrix with diagonal terms measuring the ISC, ρ_d, of the corresponding projections of the data, y_d =  v_d^⊤ x. Note that these additional projections of the data are uncorrelated, as desired, because their correlation matrix, R^y_W  — defined as in Eq. (<ref>) but for y — is diagonal:


    R^y_W =  V^⊤ R_W  V =  Λ_W .


The fact that Λ_W is a diagonal matrix is a well-know property of the generalized eigenvalue equation (<ref>) (see Appendix <ref>).

To summarize, the principal eigenvector of the between-subject over the within-subject covariance matrices provides the projection of the data that maximizes ISC. Subsequent smaller eigenvectors provide additional projections of the data that are uncorrelated with one another. The corresponding eigenvalues are equal to the ISC that each projection achieves. It follows that the last projection captures the direction with the smallest ISC. Note that the corresponding eigenvectors are uniquely defined, provided that the ISC are distinct (see Appendix <ref>). 

In the illustrative example of Figure <ref> with D=2, the eigenvalues were ρ_1=1 and ρ_2=0.04, for the red and blue directions respectively (see Figure <ref> D and <ref> F). Evidently the shared dimension where the two subjects vary in unison (red) has been identified. A concrete application of CorrCA to the problem of identifying brain activity that is shared between subjects is described in section <ref>. A thorough evaluation of the conditions under which the algorithm identifies the correct projection vectors v_d will be presented in  section <ref> (Fig. <ref>). The main finding is that the shared waveforms (signal) need to project in the same direction across subjects. Performance degrades if there is additive noise, or equivalently, if the directions of the non-reproducible waveforms (noise) differ between subjects.  

Given a set of extracted correlated components, the question arises as to how many of these components exhibit statistically significant correlations. For zero mean, independent test data we can use the F-statistic to determine statistical significance (see section <ref>, Eq. <ref>). Statistical testing is described in Appendix <ref>, where we also treat the case of non-independent data using shuffle statistics. The methods for statistical testing will be validated with numerical simulations in section <ref>.  



 §.§ Fast computation of inter-subject correlation



In Eq. (<ref>), the term k=l is excluded from the sum over k. Note that the excluded term is precisely R_W, so that adding the two covariances together completes a sum with all pairs, including k=l. We will refer to this, therefore, as the total covariance:


    R_T    =    R_B +  R_W 
    
       =   ∑_i=1^T ∑_l=1^N ∑_k=1^N 
    ( x_i^l - x̅_*^l)( x_i^k - x̅_*^k) ^⊤ .


This relationship is useful because R_T can be simplified to: 


    R_T = N^2 ∑_i=1^T  
    (x̅_i^* - x̅_*^*)(x̅_i^* - x̅_*^*) ^⊤ ,


where x̅_i^* = 1/N∑_l=1^N  x_i^l and x̅_*^* = 1/T∑_i=1^T x̅_i^* are the mean across subjects and the grand mean, respectively.  Thus, surprisingly, to compute inter-subject correlation, one never has to explicitly correlate pairs of subjects, because R_B can be computed from R_T and R_W, and neither one requires a sum over pairs of subjects. With an order of N operations, this makes the computation more efficient than the direct implementation of Eqs. (<ref>)-(<ref>), which requires an order of N^2 operations. The direct implementation was used in all previous work <cit.>, but becomes problematic when one uses populations sizes reaching hundreds of subjects or more <cit.>. 



 §.§ Relationship to Linear Discriminant Analysis


 
Perhaps the naming and the additivity of the covariance matrices introduced here — Eqs. (<ref>), (<ref>), (<ref>) — reminds some readers of the between-class and within-class scatter matrices used in Linear Discriminant Analysis (LDA). In the current notation these are:


    S_T    =   ∑_i=1^T ∑_l=1^N  
    ( x_i^l - x̅_*^*)( x_i^l - x̅_*^*) ^⊤
     S_B    =   ∑_i=1^T   N
    (x̅_i^* - x̅_*^*)(x̅_i^* - x̅_*^*) ^⊤
     S_W    =   ∑_i=1^T ∑_l=1^N  
    ( x_i^l - x̅_i^*)( x_i^l - x̅_i^*) ^⊤ .


Here, index i enumerates the different classes and index l enumerates the exemplars in each class (see Table <ref> in Appendix <ref> for an overview on the use of indices in different algorithms). Therefore x̅_i^* =  1/N∑_l=1^N x̅_i^l and x̅_*^* = 1/T∑_i=1^T x̅_i^* are the class mean and grand mean vectors, respectively. Scatter matrices also satisfy an additivity rule: S_T= S_B + S_W <cit.>. The goal of LDA is to maximize the separation between classes, which is defined as the ratio of between- to within-class variance: 


    S =  v^⊤ S_B  v/ v^⊤ S_W  v .


As before, the maximal separation S is obtained with the eigenvector of S_B^-1 S_W corresponding to the largest eigenvalue S. As before, we can extract additional projections that achieve separation of classes. If we require that these projections are uncorrelated, then they are given by the additional eigenvectors of S_B^-1 S_W. This concept was first proposed by <cit.> and generalizes the case of two classes introduced by <cit.>. We refer to this as classical LDA, and note that it maximizes the sum of variance ratios S computed separately for each projection (see  <ref>). An alternative approach is needed if instead the goal is to maximize the ratio of summed variances <cit.>, or if one prefers orthogonal projection <cit.>. 

This optimality criterion and resulting eigenvalue problem look strikingly similar to that of CorrCA in (<ref>). Despite the similar naming, however, the scatter matrices differ from the between-subject and within-subject covariance matrices, and it is not immediately obvious how the two optimality criteria relate to one another. There is, however, a close relationship between the two (see Appendix <ref>): 


    N  S_B    =    R_W +   R_B , 
    
    N  S_W    =    (N-1)  R_W -  R _B   + N T  S_M ,


where S_M is the covariance of x̅^l_* across subjects:


    S_M = 
     ∑_l=1^N (x̅_*^l - x̅_*^*)
                  (x̅_*^l - x̅_*^*)^⊤ .


When these sample-mean vectors are equal across subjects then the scatter matrices are a linear combination of the between- and within-subject covariance matrices. Note that equal sample-means, x̅_*^l=x̅_*^k, does not imply equal class-means, x̅_i^*=x̅_j^*. Zero-mean signals will naturally satisfy this condition and will typically still have different class means. Under this assumption, we can also write a simple functional relationship between class separation S and inter-subject correlation ρ:


    S = r_W+r_B/(N-1) r_W - r_B = ρ+(N-1)^-1/1-ρ.


This relationship will be useful in section <ref> where we establish a link between ISC and the classic F-statistic (Eq. <ref>), which yields a parametric test for statistical significance of ISC.  Note that S is monotonically increasing with ρ (for ρ≤ 1), because the slope of this relationship is strictly positive: 

    ∂ S/∂ρ = N/N-11/(1-ρ)^2 > 0.


This means that maximizing class separation S in (<ref>) also results in maximal ISC ρ, provided that we equate the means. Thus, finding vectors v that maximize S of Eq. (<ref>) gives the same set of solutions as maximizing ρ in Eq. (<ref>). 

To see this, note that the goal of maximizing a ratio of two matrices, such as in Eq. (<ref>), can also be achieved by joint-diagonalization of the two matrices (see Appendix <ref>):


    V^⊤ R_B  V   = Λ_B , 
     V^⊤ R_W  V   = Λ_W ,


where Λ_B and Λ_W are diagonal matrices. Note that if V jointly diagonalizes two matrices, it also diagonalizes any linear combinations of the two <cit.>. Because S_B and S_W are linear combination of R_B and R_W, then the V that satisfies (<ref>) also diagonalizes S_B and S_W.  Thus, this V captures the eigenvectors to R_W^-1 R_B as well as S_W^-1 S_B, and satisfies in both cases the constraint of mutually uncorrelated components. Because the  eigenvalues of the two problems are monotonically related, the eigenvectors will be sorted in the same order for both problems. Thus, the solutions to both optimization problems are the same, including the constraint of uncorrelated components. To summarize, CorrCA and classical LDA yield the same result, provided that the sample mean vectors are equal across subjects.

The equivalence between CorrCA and LDA is perhaps surprising. Note that LDA attempts to separate classes so that exemplars of one class do not overlap with the exemplars of other classes. In the illustrative example of Figure <ref>, time samples i take on the role of classes. There are 20 classes, each with two exemplars (connected with gray lines). Classes are well separated in the direction of component y_1 (red arrow), whereas they are overlapping in the original dimensions x_1 and x_2.  Note that the signals are zero-mean across samples but have differing class-means. In this toy example, there is no variation between the two exemplars in each class (vertical lines in panel F indicate that the y_1 values are identical within a class). Thus, CorrCA has found a dimension of the data (y_1) in which the classes are perfectly separated. Zero within-class variance leads to infinite separation (Eq. <ref>) and this in turn to unit correlation  (Eq. <ref>).  More generally, we have demonstrated with Eq. (<ref>) that increasing separation between samples is equivalent to increasing correlation across samples. 
To understand this, note that separation is high when multiple repeats have similar values, yet these values have to be different for different samples. This means that multiple repeats have to vary together across samples. In other words, the variations across samples have to be highly correlated across repeats.

Notice that the equivalence between CorrCA and LDA requires the same number of exemplars in each class – something that is not usually the case in conventional LDA. In fact, there is a one-to-one link between exemplars across classes (gray lines in panels C and F of Figure <ref>).
  









 §.§ Maximum reliability and the F-statistic



In the present context, the scatter matrices S_B and S_W are interesting for another, perhaps more important reason. Consider the case where l represents repeated measurements of raters, or repeated measures on the same subjects, so that ρ now measures inter-rater correlation, or inter-repeat correlation (IRC). According to Eq. (<ref>), S_B measures  (except for a scaling factor) the sample-covariance of x̅_i^*, which is the mean across repeats. On the other hand, according to Eq.  (<ref>), S_W measures the sample-covariance averaged over repeats. When projected onto y with v they define  the variance of the mean across repeats, σ_y̅^2, and  average variance around these means, σ̅_y^2, respectively.  The ratio of the two variances is a sensible definition for signal-to-noise ratio (SNR):


    S = σ_y̅^2/σ̅_y^2 .
 

We take this SNR also as a metric of repeat-reliability. In this view, the results of the previous section show that maximizing correlation between repeats is equivalent to maximizing repeat-reliability. In particular, Eq. (<ref>) provides the relationship between SNR and ISC. An application of CorrCA to the problem of identifying ratings that are reliably reproduced between different raters is described in Sections <ref> and <ref>. In Section <ref> we describe results on the problem of identifying components of brain activity that is reliably reproduced across repeated trials. In all instances, the goal is to maximize ρ, which is equivalent to maximizing S.

A component decomposition method that specifically maximizes SNR
was previously described in <cit.>. This method can be formulated as a joint diagonalization (JD) problem similar to CorrCA or LDA. When SNR is defined as in Eq. (<ref>), then the JD approach diagonalizes matrices S_T and S_B in the present notation. By the same argument as in section <ref>, it is clear that this JD approach provides the same solution as LDA and CorrCA. The optimization criteria for the three techniques are summarized in Table <ref>.



We note that S is a ratio of variances, which permits a direct link to the F-statistic, as used in the analysis of variance:


    F = T(N-1)/T-1 S = T(N-1) ρ+T/(T-1)(1-ρ) .
  

For normally and independently distributed samples this quantity follows the F-distribution with degrees of freedom d_1=T(N-1) and d_2=T-1  <cit.>.  We can therefore use the F-distribution to assess statistical significance for ρ. This gives us a strict parametric test of significance for each CorrCA component (on unseen data and provided the means are equalized, i.e. x̅^l_*=x̅^*_*). We will leverage this observation when we perform tests for statistical significance of ρ evaluated on independently and identically distributed test data. For non-independent samples, which are more typical of the temporal signals of interest, we will rely instead on shuffle-statistics, which can also be applied directly to training data (see Appendix <ref>). Both methods will be evaluated for accuracy in section <ref>. 
 


 §.§ Forward model



An import aspect of multivariate models is parameter interpretation. CorrCA is an example of a “backward model” in the sense that the observed data is transformed into components that capture the source of covariation. It shares this property with a host of other methods such as PCA and ICA. An alternative strategy is “forward modeling” such as Factor Analysis or General Linear Models, which transforms known variables to explain the observed data <cit.>. Such models often capture the physical processes underlying data generation. For instance, in the case of electromagnetic brain signals, a forward model refers to the “image” that a current source in the brain generates at the scalp sensors.  

For CorrCA, the backward model is given by the matrix V which maps observations x to components y= V^⊤ x. A corresponding forward model can be defined as the projection A that best recovers measurements x from the components y <cit.>: 


    x̂ =  A y.


The least-squares estimate of this forward model projection is


    Â= R_W  V( V^⊤ R_W  V)^-1.


Note that the columns a_k of this matrix capture how correlated a putative source y_k is with the observed sensor recordings x. This approach of recovering the forward model from the backward model is explained in more detail in <cit.>, along with a discussion on how forward and backward model parameters need to be interpreted. In the illustrative example of Figure <ref>, the forward model is shown in panel C as red and blue arrows for first and second component respectively. Evidently, the directions recovered by CorrCA match the underlying source orientations described in that example. 

The utility of the forward model will be demonstrated for the case of brain activity, where we are interested in how the activity of a recovered neural source manifests at the sensor level (Figure <ref> and Figure <ref>). These forward models can be used directly to perform source localization to identify the spatial origin of the corresponding current sources <cit.>. In the case of questionnaire ratings (section <ref>) we will interpret the meaning of a component by inspecting the forward model, as it measures the correlation of a component with answers of raters on different topics.  
For the case of ratings, we may also be interested in which items lead to a reliable aggregate score (sections <ref>). In these cases, an analysis of backward model weights may provide valuable insight on what data features are most important to construct a reliable aggregate measure.

Note that if the original data is uncorrelated, then R_W is diagonal and the projections V are orthogonal. In that case the forward and backward models are identical, except of an overall scale for each component. In that case, we can directly inspect the backward model, as is customarily done in PCA, which has orthogonal projections and one inspects the component “loadings”. In the example of Figure <ref>D the covariance R_W is approximately diagonal and one can inspect the results with either the forward of backward model.




 §.§ Relationship to multi-set Canonical Correlation Analysis




Thus far we have assumed that all repetitions should receive a common projection vector. This is sensible in the case of multiple trials if we expect that the sources of interest behave similarly across repetitions. However, in the case of diverse subjects or raters, it may be more appropriate to assume that each repetition should be uniquely weighted. For instance, brain activity recorded at a given location may differ between subjects due to anatomical or functional differences. This brings us to the case treated by Canonical Correlation Analysis (CCA) for M=2, or more generally, multi-set CCA (MCCA) for M>2. There are multiple implementations of MCCA <cit.>. We will derive one of these instantiations as a maximization of ISC as defined in Eq. (<ref>).  Assign to every subject l a unique projection vector v^l:


    y^l_i =  v^l⊤ x^l_i .
 

By inserting this into definition (<ref>), it follows readily that


    ρ =  
    ∑_l=1^N ∑_k=1, k ≠ l^N  v^l⊤ R^lk v^k
    /
    (N-1) ∑_l=1^N  v^l⊤ R^ll v^l
     ,


where R^lk is the cross-covariance matrix between x_i^l and x_i^k:


    R^lk 
    =
    ∑_i=1^T 
    ( x_i^l - x̅_*^l)( x_i^k - x̅_*^k) ^⊤.



The projection vectors v^l that maximize ρ can be found as the solution of  ∂ρ / ∂ v^l⊤ = 0, which yields for each l the following equation: 


    1/N-1∑_k=1, k ≠ l^N  R^lk v^k
    =
    ρ  R^ll v^l   .


This set of equations for all v^l can be written as a single equation via concatenation to a single ND-dimensional column vector v = [  v^1⊤… v^N ⊤]^⊤:


    R v
    =
    λ  D v ,


where, λ=(N-1)ρ+1, R is a matrix combining all R^lk, and where the covariance matrices R^ll  are arranged in a block-diagonal matrix D:


    R = 
    [
    [ R^11 R^12    … R^1N; R^21 R^22    … R^2N;    ⋮    ⋮    ⋱    ⋮; R^N1 R^N2    ⋯ R^NN;      ]], 
     D = 
    [
    [ R^11    0    …    0;    0 R^22    …    0;    ⋮    ⋮    ⋱    ⋮;    0    0    ⋯ R^NN;      ]].


If D is invertible, the stationary points of the ISC are now the eigenvectors of D^-1 R, or more generally, the generalized eigenvectors of Eq. (<ref>). We can arrange all such eigenvectors as the columns of a single matrix V = [ v_1 … v_n … v_D]. The eigenvector with the largest eigenvalue λ maximizes the ISC because ρ increases linearly with λ. The subsequent eigenvectors maximize ISC subject to the constraint that the component signals are uncorrelated from the previous components (see Appendix <ref>).  In other words, the solution V diagonalizes D:


    V^T  D V = Λ.
  

Formulating the MCCA problem as an eigenvalue problem can be derived in multiple ways  <cit.>. In <cit.> it is presented as the maximization of the eigenvalue of the correlation matrix between subjects (MAXVAR criterion), but it can also be derived by maximizing the sum of correlations between subjects (SUMCORR criterion) <cit.>. Here we have derived MCCA as a maximization of ISC. In this sense CorrCA can be seen as a version of MCCA constrained to have identical  projection vectors for all N data sets (repeats/subjects). MCCA has ND parameters per component, whereas CorrCA has only D. When there is an abundance of data (T ≫ ND) and there is reason to believe that the N datasets should be weighted individually, then MCCA may be appropriate <cit.>. However, when N is large compared to the data size  (N>T/D) CorrCA may be preferable due to the reduced number of free parameters. This will be demonstrated in section <ref>. 

When there are only two data sets (M=2) then the present MCCA algorithm reduces to conventional CCA <cit.>. In the illustrative example of Figure <ref>, CorrCA by definition gives the same result for both subjects (red and blue arrows), but we find that the directions derived with CCA differ between subjects (not shown). In section <ref>, we will apply CCA to ratings data from two distinct raters (pairs of parent and child), in which case distinct projections may be warranted.




§ APPLICATIONS


In the following sections we demonstrate existing and new applications of CorrCA, while also examining issues related to statistical significance testing and regularization of the model. We also compare CorrCA with MCCA/CCA and kernel-CorrCA. 



 §.§ Extracting brain activity that correlates between subjects



We have previously used CorrCA to analyze brain activity recorded from multiple subjects  <cit.>.  Briefly, subjects in these experiments watched movies or listened to narrated stories while their neural responses were recorded with multiple sensors (electroencephalography (EEG) — electrical potentials measured at multiple locations on the scalp surface).  CorrCA was then used to extract components that were most reliable across subjects. These experiments showed that brain signals are significantly correlated between subjects on a rapid time-scale of less than a second (signals were high-pass filtered at 0.5 Hz). While the ISC values are small (ρ≈ 0.05), these values are statistically significant and very reliable across sessions <cit.> even in realistic environments <cit.>. Interestingly, the level of correlation between subjects as measured by ISC is indicative of how attentive subjects are towards the stimulus <cit.> and, therefore, predictive of how well individuals remember the content of the videos <cit.>.  ISC of the EEG also predicts the retention of TV and online audiences <cit.>. 






Here we reanalyzed the data of <cit.> to compare CorrCA with MCCA (see section <ref>) on their ability to learn reliable projections of the data, while also using the opportunity to demonstrate forward model computation (section <ref>) and statistical testing (Appendix <ref>).   The data set consisted of the EEG from N=18 subjects as their brain activity was captured with D=64 scalp electrodes across T=50344 time samples corresponding to the 197-second video stimulus.  We first compared CorrCA with regularized MCCA in terms of the number of reliable dimensions extracted by each, computed here using a non-parametric test using circular shuffle statistics (see Appendix <ref>).   The ISC values found on the entire data set are shown for the first 5 components in Figure <ref> A , where blue bars correspond to CorrCA and green to MCCA.  Moreover, darker shading indicates that the ISC of the component was found to be statistically significant: 3 significant dimensions were found by CorrCA, while only 1 was found by MCCA, despite the (training set) ISC values being higher for the latter.   
Indeed, when evaluating ISC on unseen test data the ISC drops substantially for MCCA as compared to the training data, and CorrCA yielded higher ISC values than MCCA for 4 of the first 5 dimensions (Figure <ref> B). 

Both of these findings are indicative of over-fitting with MCCA, despite the use of regularization. We used the Truncated Singular Value Decomposition (TSVD, similar to what is described in Appendix <ref>) and selected a subspace of K=12 dimensions as this gave the largest test-set ISC (summed over the first 3 components). Note that MCCA fits a model for each subject. In this example it uses 216 parameters per component (N=18 subjects, K=12 dimensions), whereas CorrCA uses D=64 parameters per component on these data. While stronger regularization of MCCA is possible, the regularization is agnostic to the shared dimensions, which explains why ISC dropped with stronger regularization on the test data. 
CorrCA can be seen as a form of regularization that addresses over-fitting, without removing dimensions that significantly correlate between subjects. 

Both CorrCA and MCCA learned spatial filters that were applied to the multi-electrode EEG data  to maximize the ISC.  The activity captured by each spatial filter (component) is best depicted by “forward models” which are illustrated row-wise for the first three correlated components in Figure <ref> C & D.  The components extracted with CorrCA are marked by a smooth topography consistent with activation of the cerebral cortex (Figure <ref> C). The components produced by MCCA are less smooth (Figure <ref>D shows subject-specific forward models for the first three components of 5 of the 18 subjects, i.e., subjects 1, 5, 9, 13, and 17).  Note that only the first component replicates, for some subjects, the bilateral occipital activation indicative of visual cortex activity seen in panel C (top). Other topographies are unlikely to be generated by a genuine neural source, and instead further support the assertion that MCCA overfit the data. 

Here and in the following examples we computed model parameters v by leaving out samples i, and measured test-set ISC on the left-out samples. The alternative could have been to leave out subjects. However, the number of entries in the between-subject correlation matrix R_B defined in Eq. (<ref>) scales quadratically with the number of subjects N, but linearly with number of samples T. Thus, we expect the accuracy of R_B and the resulting eigenvectors v in Eq. (<ref>) to  decrease faster when removing subjects instead of samples.


  




 §.§ Maximizing inter-rater reliability



The theory laid out in section <ref> shows that one can use CorrCA to identify components with maximal inter-rater reliability. In this following example, clinicians assessed the motor skill of children by measuring performance on various standard tasks (Zurich Neuromotor Assessment). The rating obtained for each task here is the time taken by the child to perform each task. This time is normalized relative to the values of a standard group of children of the same age <cit.>. Although this is an objective measure, there is nonetheless variability due to variable performance of the child when repeating the task, and also subjective variability in the observer operating the stop-watch. In this specific data set, we have ratings from T=30 children on D=12 tasks from N=2 raters. The 12 measures are mostly positively correlated across the children (Figure <ref> A, see panel G for descriptive task labels). Reliability between raters is measured here as the IRC, which differs across tasks (Figure <ref> B). 




For many assessment instruments, individual ratings are summed up to obtain a total score. It is also natural to compute sub-scores by combining ratings on related items. In these sums, ratings effectively obtain equal weights, which is a sensible but nonetheless arbitrary choice. CorrCA can provide weightings that will give the most reliable response across raters, as reproducibility of results is a desirable goal.  In this specific example, the goal is to find the weightings of the 12 motor tasks that provide reliable aggregate (sub-) scores, and which capture independent aspects of these ratings (independent to second order). The resulting weighting vectors V are shown in Figure <ref> G. Each column corresponds to a different component and has been normalized to have unit norm. 
Components are sorted in decreasing order of IRC (Figure <ref> E) and therefore decreasing order of inter-rater reliability (see section <ref>, Eq. <ref>).  Of these components, the first five reach statistical significance at an α-level of 0.05 (highlighted in Figure <ref> E). Inspecting the “backward model” V it is clear that first component relies heavily on a single task (task 2, which has the highest IRC) but has small positive and negative contributions from other tasks. Components 2 and 3 are mostly a contrast between tasks 3 and 8, and 2 and 4, respectively. The detailed results can be obtained by running the code accompanying this article. In total, the method extracts different aggregate scores of high reliability, relying in some instances on a small subset of tasks. This has the potential of simplifying the assessment of motor performance while increasing its reliability.  An important aspect of these various components is that, by design, they are uncorrelated from one another (uncorrelated across the 30 children; Figure <ref> D), whereas the original task ratings all measure similar aspects (they are predominantly positively correlated; Figure <ref> A).

In this data set, children were also evaluated a second time by the same rater in order to assess test-retest reliability (Figure <ref> C). Comparing panels B and C, it is evident that the test-retest reliability is somewhat higher than the inter-rater reliability (i.e., the evaluation does indeed depend to some degree on the subjective judgment of the rater). The components extracted with the highest IRC also exhibit the strongest test-retest correlations (at least for the first few components). This confirms that these aggregate measures generalize to data that was not used in the CorrCA optimization, and that minimizing subjective variability of the observer also minimizes test-retest variability. 





 §.§ Identifying agreement between child and parent



We also provide an additional example of CorrCA on subjective ratings data related to the style of parenting. In this example, the goal is to identify questions on which parent and child agree. Specifically, we analyzed data collected with the Alabama Parenting Questionnaire <cit.> as part of the Healthy Brain Network initiative <cit.>. In the APQ, parents and children answered D=42 questions pertaining to their relationship by providing a numerical value from 1 to 5 for each question. From this, the APQ averages various questions to provide aggregate scores related to different aspects: “involvement score”, “positive parenting score”, “poor supervision score”, “inconsistent discipline score”, and “corporal punishment score”.  At the time of this analysis, complete data were available for T=616 children/parent pairs (ages 6-17, 40% female), i.e. N=2 raters. The first observation from the data is that the answers from children and parents are only weakly correlated (ρ=0.17 ± 0.1, Figure <ref> A). The aggregate scores fare even worse (ρ=0.064 ± 0.13). A two-sample t-test suggests that these aggregate scores are less reliable than individual questions (t(47)=2.4, p=0.022) including several scores with negative correlations. Evidently, children do not much agree with their parents when it comes to parenting issues, despite answering nearly identical questions on the topic.  

We applied CorrCA to find a combination of questions on which parents and children agree. To this end, we divided families at random into a training and a test set of equal size (T=308) and computed the projections V that maximized IRC on the training set. We then evaluated the resulting projections in terms of the correlation between child and parent for both train and test sets (Figure <ref> B & C). By design, IRC on the training set is largest for the first component. On the test set, the first component also has a relatively large IRC, whereas the other components do not differ substantially in IRC from the original set of questions (IRC hovers around 0.2 in the original ratings, X, and the aggregate component scores, Y). The lower IRC values on the test set as compared to the training set indicate the presence of over-fitting, although the discrepancy has been reduced here using shrinkage regularization (using γ=0.5; see Appendix <ref>). To determine how many significantly correlated dimensions are actually present in these data, we used circular shuffle statistics (see Appendix <ref>) on all the data (T=616) and found K=8 components with p<0.05. On the test data alone, computation of the F-statistic revealed K=8 significant components (Figure <ref> B & C, dark shaded).



The component with the highest IRC has a combination of questions contributing to reasonably reliable aggregate scores (ρ=0.61 on test set, Figure <ref> G).  To interpret this component we inspect its forward model weights (see section <ref>.)
The following six questions correlate most strongly to the first component  (weightings of V are given in parenthesis): 

Question 32: “You are at home without an adult being with you” (-0.38); Question 11: “Your parent helps with your homework” (+0.35); 
Question 19: “Your child goes out with a set time to be home” (-0.33); 
Question 21: “You go out after dark without an adult with you” (-0.33); Question 41: “You use time-out as punishment” (+0.33).
Question 6:  “You fail to leave a note or let your parents know where you are going” (-0.28). 
The sign of the forward weights indicate here whether the answer to that question correlate positively or negatively with the aggregate measure <cit.>. Therefore, this measure seems to capture the independence of the child versus parental involvement. The positive correlation of Question 41 with this independence-vs-involvement dimension is interesting as it suggests that this form of punishment is associated with a child's independence. In total, while parents and children largely disagreed, there is moderate agreement on how independent their lives are. Given this result, we suspected that this metric increases with the child's age. Indeed, there is a strong correlation of this component with the child's age across the cohort (Pearson's r=0.7 and r=0.67 for parent and child scores respectively, p<10^-80, T=616; with no significant difference with the child's sex).  We did not analyze subsequent components due to the proximity of residual eigenvalues (Figure <ref>, panel C) and noise level may not permit accurate identification (see section <ref>). 

Given that parents and children may genuinely differ in their judgments, it would make sense to allow for different weighting to the questions.  Therefore, we also applied MCCA to this data as described in section <ref>. We found only K=4 significant components at p<0.05 (as opposed to K=8 for CorrCA; Figure <ref>E, dark shaded). ISC summed over these first four components was numerically smaller for MCCA vs CorrCA (0.69 vs 1.10, values obtained on the test data without regularization; Figure <ref> F). Therefore, it appears that on this data the sample size of T=616 was not sufficiently large to support a doubling in the number of free parameters (from 42 for CorrCA to 84 for MCCA). 

We have also tested the non-linear kernel-CorrCA method (see Appendix <ref>) on this ratings dataset (Figure <ref> E). There was a small gain in terms of inter-rater correlation or inter-trial correlation respectively in the test data (Figure <ref> F). However, the projections of the strongest components were not substantially different from the linear projections, suggesting that there were no strong non-linear relationships in the data. We find a similar result with kernel-CorrCA for the the data presented in the next section <ref> (results not shown).   



 §.§ Identifying stimulus-related components in EEG with high SNR



PCA is often used to extract a low-dimensional representation of high-dimensional data. To this end, PCA selects orthogonal projections that capture most of the variance of the projected samples. When multiple repetitions of the data are available, CorrCA can be used for the same purpose of reducing dimensionality of the data, but with the objective of capturing dimensions that are reliably reproduced across repetitions. Here we apply this to a conventional scenario in neuroimaging: neural activity is collected over time with multiple sensors, while the same experimental stimulus is repeated in multiple trials. We propose to use CorrCA to identify spatio-temporal patterns of activity that are well preserved across trials, similar to what was proposed in <cit.>. Ideally, the corresponding time-courses of multiple components are uncorrelated in time so that they capture different activity. CorrCA guarantees this by design as it diagonalizes R_W. CorrCA is therefore similar to other source separation methods such as Independent Components Analysis <cit.>, Denoising Source Separation <cit.>, or Joint-Decorrelation  <cit.>. In fact, CorrCA can be viewed as a form of blind-source separation, as no labels are required for the analysis.


Here we analyze EEG data that were collected during an Eriksen Flanker Task <cit.>. Subjects were presented with a collection of arrows pointing left or right. They were asked to decide as quickly as possible between two alternative choices (the direction of the center arrow) by pushing corresponding buttons with either their left or right hand. We analyze data from a single participant who responded correctly in 550 trials, and made an error in 46 trials (Figure <ref>). This error, which participants notice as soon as they make their incorrect response, leads to a negative EEG potential 0-100 ms after the button press (this period is indicated between black and red vertical lines in panel B & C for electrode #7, see electrode location in panel H). This phenomenon is known as the “error-related negativity” <cit.>.   
Button pushes also elicit the “readiness potential” — a gradually increasing potential preceding the button push. Aside of the readiness potential and error-related negativity there are a series of other potential deflections following the button push. In this example, we have T=176 time samples, D=64 electrodes, and N=46 trials. Time is aligned across trials so that the button press is at time t=0 ms. The traditional approach to analyzing such data is to average the electric potentials over trials, yielding what are called event-related potentials (ERPs, panel B). The reliability of these trial-averages is relatively low (panels B& D show mean and standard-error across trials, and panels C & E show single trial data). We selected  electrodes #7 and #22 for display as these capture distinct time courses and locations identified in prior literature. The SNR (as defined in section <ref>) is maximal for the same cluster of frontal electrodes (#12, #17, #44, #43) and the mean SNR across all electrodes is S√(N) = 0.72 (Figure <ref> A). We are multiplying SNR with √(N) because, typically, the goal is to determine if the mean value differs significantly from zero, and thus standard-error of the mean (SEM) is the relevant quantity. 




In our earlier work we specifically looked for a projection that distinguishes the EEG evoked response between error and correct responses <cit.>. Here the goal is to identify reliable error responses. When we apply CorrCA to these data, we obtain a series of components with high inter-trial correlation (Figure <ref> F). CorrCA components were extracted using shrinkage regularization (γ=0.4, see Appendix <ref>) on a training set of N=45 trials (leaving 1 trial out for testing). Performance was evaluated on the leave-one-out test set. The first three components (highlighted in dark blue) are significant at the α=0.05 level according to a non-parametric test based on circularly shifted surrogate data (see Appendix <ref>). The first two of these components have an SNR of S√(N) = 3.3 and 1.6. Evidently, component 1 is very reliable (panel G shows that the mean is several SEM away from the baseline). The deflections are now obvious even in single trials without the need to average (panel I). The spatial distribution of this component activity (panel H) is indicative of the well-known ERN. Component 2 (panel J) is indicative of the readiness potential <cit.>. It rises leading up to the button push at time t=0 ms, and has a spatial distribution over parietal and motor areas (panel K).

Note that CorrCA has reduced the D=64 dimensional data to just a few components that capture the reliable portion of the data (high SNR; panel F).  By compactly representing neural activity in a small number of components, it becomes easier to probe for changes in activity due to manipulation of experimental variables, whose effects on the original D=64 dimensional data may be obscured due to the lower SNR and the multiple comparison problem.  

Like most data decomposition algorithms, CorrCA does not require that the underlying data have any time structure. For time series data, this implies that the correlated signal of interest is observed with the same delay in each channel and repetition. If this assumption is violated, drops in performance will occur.



 §.§ Effect of regularization on the reliability of EEG responses



A common numerical difficulty with generalized eigenvalue problems is the required matrix inverse. In the case of CorrCA, it is necessary to invert R_W in equation (<ref>). When some of the eigenvalues are small (or even zero), then random fluctuations due to noise will dominate the inverse. It is thus important to regularize the estimate of R_W. In the past, we have used the Truncated Singular Value Decomposition (TSVD) <cit.> and shrinkage <cit.> for this purpose <cit.>.  These methods are described in detail in Appendix <ref>. 

To evaluate the effects of these regularization methods on the results of CorrCA we use data from our study on neural responses to advertisements <cit.>.  For both the TSVD and shrinkage approaches, we spanned the space of regularization strength (a single parameter, K and γ, in each case) and measured ISC.  The training set was the EEG data from 12 subjects during the first viewing of the ad, while the EEG data collected during the second viewing served as the test set.  

The sum of the first 3 ISC values for both regularization methods followed a similar trajectory with increasing regularization strength (decreasing K and increasing γ; Figure <ref> B & C). On the training data, ISC drops monotonically, whereas on the test data a maximum is achieved at intermediate regularization values. On these data, shrinkage regularization was most effective at γ=0.4 and TSVD when retaining K=20 out of a total of D=64 dimensions in the within-subject covariance R_W.  The forward models of the first three correlated components for each regularization technique are shown in Figure <ref> A.  Notice that regularization tends to smooth out the topographies of the second and third components compared to no regularization.    





§ VALIDATION ON SIMULATED DATA



We evaluated the conditions under which CorrCA identifies the true source directions. We also determined the performance of CorrCA at identifying components with high ISC, and compared this against a simple benchmark based on PCA of the subject-average data. Additionally, we compared the accuracy of the statistical tests described in Appendix <ref> at estimating the correct number of components. To this end, we generated artificial data with known properties. The data are generated assuming a known number of correlated components (signals shared across subjects) and additive spatially correlated noise not shared across subjects. The following factors were varied parametrically: the number of subjects (raters, repeats), N, the number of samples per subject, T, the number of dimensions, D, the number of correlated components, K, and the signal-to-noise ratio (SNR) as defined in the space of multivariate measurements. Each simulated correlated component was identically reproduced in each subject, so that the ISC of all components was equal to one. We further varied the temporal dependency structure of the signal and noise components (IID or pink noise), the distribution of the signal and noise components (Gaussian or χ^2 distributed), and the heterogeneity of the signal and noise subspaces (either the same or different for all subjects). Signal components were generated to be correlated across subjects, but uncorrelated with each other. Noise components were uncorrelated with each other as well as across subjects. The number of noise components was fixed to D. The following default parameters were used unless otherwise noted: N=5, T=200, D=30, K=10, SNR = 0 dB, no regularization, IID Gaussian signal and noise components, ISC = 1 for all K components. Details on data generation and evaluation methods can be found in Appendix <ref>. 




 §.§ Accuracy as a function of Signal-to-Noise Ratio



Figure <ref> depicts CorrCA performance as a function of measurement SNR for fixed T=200, D=30, N=5, K=10. SNR values range from -40 dB to 40 dB. It is apparent that average training and test ISC reach the optimal value of 1 for very high SNR (panel A). For very low SNR ranges, training ISC significantly exceeds test ISC, indicating over-fitting. The identification of the correlated subspace converges to the optimal value for high SNR (panel B). However, individual components are not perfectly identifiable even for high SNR because they share the same ISC values (see Appendix <ref> for a discussion). The identifiability of individual components becomes possible if the ISC of the simulated signal components are adjusted to linearly decrease from 1 down to 1/K, and if the measurement SNR is high (data not shown). In practice, one should keep in mind that the uniqueness of the ISC (and thereby, their individual identifiability) depends on the amount of noise that is collinear to each original correlated dimension. In the IID case, all estimates of K converge to the true number of 10 for high SNR (panel C). Here, the parametric test and the surrogate data approach using random circular shifts perform similarly, while the surrogate data approach using phase scrambling performs considerably worse at low SNR, presumably due to the absence of a well-defined notion of frequency and phase in the sense of the Fourier transform. In the case of dependent (pink noise) samples, both surrogate data approaches perform identically, converging to the true value of K for high SNR (panel D). In contrast, the parametric test, which incorrectly assumes IID data, overestimates the number of correlated components across the entire SNR range, as expected. 

Note that the significance test is applied here to multiple components (D=30). Under appropriate conditions, all three methods converge to the correct number of dimensions in these simulate data (K=10). This indicates that all three methods (described in Appendix <ref>) properly account for these multiple comparisons. 






 §.§ Accuracy of estimating dimensionality of correlated components


The performance of CorrCA as a function of the true number of correlated components, K, is depicted in Figure <ref>. Here, T=200, D=30, N=5, SNR = 0 dB, where K was varied from 3 to 28. As the overlap between signal and noise subspaces increases with K, the average ISC of the reconstructed correlated components drops (panel A). A similar decline is seen in the reconstruction of simulated correlated components and their forward models (data not shown). As a result, the percentage of correlated components that can be recovered (i.e., reach statistical significance) also decreases with increasing K  (panels B & C).

We also compared the ISC values obtained with CorrCA against a simple benchmark, namely, PCA of the subject-averaged data. While this simple method captures the correct directions for high SNR case, CorrCA significantly outperforms at medium and low SNR values (Figure <ref> A). The same is true for simulations with varying number of underlying dimensions (Figure <ref> A), where we see that CorrCA outperforms this simple method in all instances tested.  




Figure <ref> shows the estimated number of correlated components as a function of the number of data points, T, and the number of repetitions (subjects/raters), N, with fixed D=30, K=10, SNR = 0 dB. In case of IID data, all estimates converge to the true value of 10 for increasing T or N. This is also true for pink noise data in combination with non-parametric statistical tests. For the parametric test, however, the number of seemingly correlated components diverges as T increases, and is overestimated regardless of N.





 §.§ Dependence on signal and noise distributions


Figure <ref> A & B shows the performance of CorrCA for different signal and noise distributions. Here, T=200, D=30, N=2, K=10, SNR = 0 dB. Signal and noise components were generated either as Gaussian or squared Gaussian (χ^2-distributed). For Gaussian distributed components, we considered an additional case, where observations x^l(t) were dichotomized (thus, Bernoulli distributed) by taking their sign. This simulation scenario resembled the typical setting in which two raters provide a binary assessment of the same issue. While non-Gaussianity of the underlying components does not substantially affect the extraction of correlated components, the same is considerably more difficult in the presence of dichotomous observations. An additional analysis was conducted for IID Gaussian data and N=5 to assess the performance of CorrCA when signal and noise covariance structures (as defined by the mixing matrices A_s^l and A_n^l, n = 1, , N) differ across subjects (Figure <ref> panel C & D).  Significant drops in the number of estimated correlated components as well as the average ISC are observed if the signal directions, the noise  directions, or both, differ across subjects.  Note that varying directions of the non-reproducible component A_n^l across subjects is equivalent to adding noise, so that the dependence on the noise directions shown in Figure <ref> (panel C & D) is a restatement of the dependence on SNR of Figure <ref>.  The result illustrates that identical directions for the  signal components across subjects are a fundamental assumption of CorrCA. If this assumption is suspected to be violated, the use of MCCA (section <ref>) may be beneficial, or approaches that can gradually move from CorrCA to MCCA <cit.>.






§ CONCLUSION


The goal of this work was to provide a formal, yet didactic and comprehensive treatment of this new component analysis technique. The analytic development resulted in a direct link to MCCA and a surprising equivalence to LDA, thus linking the new methods with classic concepts of multivariate analysis. We also found an efficient scheme for computing ISC/IRC which surprisingly does not require computation of pairwise correlations. We also identified  the F-statistic as an exact parametric test for statistical significance, which is valid in the case of normal and identically distributed samples. For the case of dependent samples we demonstrated and validated circular-shift shuffle statistics as an efficient non-parametric test of significance. 
We also presented an extension to non-linear mappings using kernels. 
Finally, perhaps the most important contribution of this work is the formalization of inter-subject correlation. Conventional metrics, such as pairwise Pearson's correlation are limited to two signals. Inter-subject correlation as defined here is applicable to more than two signals. ISC differs from conventional definitions of correlation (see section refapdx:ICC), yet seamlessly integrates in well established mathematics related to linear discriminants, canonical correlation, and the F-statistic. 

We make all code and data available at <http://parralab.org/corrca>,  and for long-term storage here <cit.>. This code is significantly faster than previous instantiations and has been extensively evaluated. All figures presented in this paper can be reproduced with this new code base and associated data. We hope that this analysis and code inspires new uses for correlated component analysis. 



§ ACKNOWLEDGMENTS


We are grateful to Valentin Rousson who pointed us to work on inter-class correlation, and Oskar Jenni of the Children University Hospital in Zurich, Switzerland who collected  the data published in <cit.>, which we reanalyzed here in Figure <ref>. Similarly, we would like to thank Mike Milhalm for pointing us to the problem of identifying parent/child agreement and Lindsay Alexander for providing access to the Healthy Brain Network data <cit.> for Figure <ref>. We also thank Samantha Cohen for providing the data from <cit.> for Figure <ref>. We thanks Samantha Cohen, Jens Madsen Michael X Cohen, and two anonymous reviewers for useful comments on an earlier version of this manuscript. We want to thank Flaviano Morone for assisting with the proof in Appendix <ref>. 








§ GENERALIZED EIGENVECTORS FOR TWO SYMMETRIC MATRICES 




 §.§ Joint diagonalization



Here we will review the well-know relationship between the joint diagonalization of two symmetric matrices A and B and the eigenvectors of B^-1 A. Consider eigenvectors arranged as columns in V with the corresponding eigenvalues in the diagonal matrix Λ: 


    A V =  B VΛ .


To find the solution V for this equation, we can replace the positive definite matrix B with its Cholesky factorization


    B= L L^⊤.


Insert this into (<ref>) and left-multiply with L^-1 to obtain:


    ( L ^-1 A L^-⊤ ) ( L^⊤ V)    =    ( L^⊤ V ) Λ , 
     M U   =     UΛ ,


with M =  L ^-1 A L^-⊤ and U =  L^⊤ V. Equation (<ref>) is a conventional eigenvalue equation with the same eigenvalues as (<ref>). It is now easy to see that the corresponding eigenvectors V =  L^-⊤ U diagonalize matrix A


    V^⊤ A V =  U L^-1 A L^-⊤ U=  U^⊤ M U = Λ .

They also diagonalize B, which we see by left-multiplying (<ref>) with V^⊤ and right-multiplying with Λ^-1. This yields:


    V^⊤ B V  =  V^⊤ A VΛ^-1 =  I .


Thus, the eigenvectors that solve (<ref>) jointly diagonalize A and B. The reverse is also true. A matrix V that satisfies:


    V^⊤ A V   = Λ_A 
     V^⊤ B V   = Λ_B ,


where Λ_A and Λ_B are diagonal matrices, also satisfy (<ref>). To see this, solve (<ref>) for V^⊤ and substitute this into (<ref>), which yields:


    Λ_B  V^-1 B^-1 A V  = Λ_A ,


and therefore 

    B^-1 A V  =  VΛ_B^-1Λ_A ,


which is the same as the eigenvalue equation (<ref>) with 


    Λ = Λ_B^-1Λ_A .




 §.§ Linear combinations



A matrix V that simultaneously diagonalizes two matrices A and B, as in (<ref>)-(<ref>), will also diagonalize any linear combination of these two matrices <cit.>, because


    V^⊤ (α A + β B)  V   =   α V^⊤ A V  + β V^⊤ B V
       =   αΛ_A + βΛ_B ,


which is obviously a diagonal matrix. Therefore, the eigenvectors of A^-1 B are also eigenvectors of (α A+ β B)^-1(γ A+ δ B), assuming vectors [α, β] and [γ, δ] are not collinear.  



 §.§ Identifiability



An additional important observation is that the eigenvectors of (<ref>) are only uniquely defined for unique eigenvalues. If an eigenvalue is degenerate, say with multiplicity K, then the corresponding eigenvectors span a K-dimensional subspace. Any vector within that subspace is a solution to the eigenvalue equation with that eigenvalue. In this space any arbitrary rotation will yield a valid solutions. Thus, when estimating components, it is possible that solutions with close-by ISC are “mixed”. Such components are not uniquely identifiable. One should also note that the sign of any eigenvector is arbitrary. If v_d is a solution, then so is - v_d. One can always arbitrarily swap the sign of v_d for any component d, which also reverses the sign of y_d^l and a_d — the dth column of the forward model A. Thus, keep in mind that in all figures the sign of individual v_d and a_d are arbitrary.  The same problem arises with most signal decomposition methods such as PCA, ICA or CCA. A common approach used in practice to “fix” the sign is to select one sensor/dimension with the strongest forward-model coefficient, and set that to a positive value (if indeed it is significantly larger than other coefficients with opposite sign).




 §.§ Optimization of sum ISC



The optimization criterion (<ref>) has the form


    J( V) = 
    ∑_d v_d^⊤ A v_d/ v_d^⊤ B v_d

and is subject to the constraints that components are uncorrelated (relative to B):

    v_d^⊤ B v_c =0, for  d  c.


It will be convenient to convert this optimality criterion into conventional Rayleigh quotients. To do so use the Cholesky decomposition (<ref>) as before and replace, v_d =  L^-⊤ w_d and A =  L M L^⊤ to yield: 

    J( W) = 
    ∑_d w_d^⊤ M w_d/ w_d^⊤ w_d

with the constraint now being a conventional orthogonality constraint:

    w_d^⊤ w_c =0, for  d  c.

One can now optimize for matrix W with the w_d as its column vectors. The solution to this constrained optimization problem is well established <cit.>. It is given by, W= U, the eigenvalue matrix of equation (<ref>). Inserting U= L V into (<ref>) gives the solution to the original problem, which is given by the generalize value equation (<ref>). 


To see that the eigenvector matrix of M maximizes J( W), first note that except for the constraint, the terms in the sum can be optimized independently, and that they all have the same form. One can therefore consider the maximum of a single term, written here simply as

    ρ =  w^⊤ M w/ w^⊤ w  ,

which is to be optimized subject to the orthogonality constraint (<ref>). From here on the argument proceeds conventionally. Let u_i, i=1 … D be the orthonormal eigenvectors of matrix M with eigenvalues λ_1 ≥λ_2 ≥…≥λ_D. We can express the solution vectors in this orthogonal basis as


    w = ∑_i=1^D α_i u_i ,
 
where α_i is the projection of w on that orthonomal basis of these eigenvectors, α_i =  u_i^T  w.  Inserting (<ref>) into (<ref>) gives 

    ρ = 
    ∑_i ∑_j α_i  u_i^⊤ M u_j α_j /∑_i ∑_j α_i  u_i^⊤ u_j α_j 
    =
    ∑_i λ_i α_i^2/∑_i  α_i^2 ≤λ_1∑_i α_i^2/∑_i  α_i^2   =
    λ_1  
    .
 
This upper bound is attained for α_1>0 and all others equal zero, and thus w_1= u_1. For the second vector we require that it is orthogonal to the first, and thus we have to constrain α_1=0. The sum in (<ref>) now excludes the first dimension and has therefore an upper bound of λ_2 which is attained for w_2= u_2. The argument continues iteratively until for the last component only, w_D= u_D, is orthogonal to all other projections. It also happens to be the minimum of ρ.  






 §.§ Uncorrelated components vs orthogonal projections



We have required here that the extracted components are uncorrelated, i.e. the projection vectors diagonalize R_W: V^⊤ R_W  V =  I. This means that the component dimensions of y_i^l are uncorrelated, but only in the averaged over all subject/repeats l. Note that for an individual subject, y_i^l may not be uncorrelated. The same is true for MCCA as derived in section <ref>. Instead of uncorrelated, we could have required that the projection vectors are orthogonal, V^⊤ V =  I, i.e. the linear transformation is limited to rotations. To find the orthogonal (rotations) that maximize ISC, one will have to resort to constrained optimization algorithms <cit.>.  A related issues arises in the context of Canonical Correlation Analysis which attempts to find a linear mapping between two datasets such that the components are uncorrelated <cit.>. Alternatively, the Procrustes problem aims to find a mapping between two data-sets with orthogonal projections <cit.>.  <cit.>. CCA and Procrustes are related, and in fact one can smoothly titrate between orthogonal and uncorrelated solutions <cit.>. In our case we favor uncorrelated components, because we want the components to capture independent information. At a minimum, this means the components have to be uncorrelated. In that sense, our method can be seen as a source separation method, similar to Independent Component Analysis <cit.>, Denoising Source Separation <cit.>, or others <cit.>.
One advantage of restricting the transformations to rotations is that euclidean distances in the new space are preserved. This means that differences found in the original sensor space are preserved in the new component space. Another advantage is that all original dimensions contribute equally to the new dimensions. In contrast, the non-orthogonal shearing operation of CorrCA or CCA may discount some dimensions and emphasize others. Both these advantages of orthonormal mappings, however, require that all new dimensions are considered. This is typically the case in the applications of hyper-alignment <cit.>. For the purposes of identifying dimensions that are reliably reproduced, however, we find in practice that many dimensions are not significantly correlated between repeats, and one will typically want to remove those. Similarly, one will want to dismiss sensors in the original data that are not reproducible across repeats while emphasizing others. Ultimately, however, the application will dictate whether one should enforce uncorrelated components or orthogonal projections. In either case, the goal can remain to maximize ISC, and the solution can be found, either through an eigenvalue question (as we presented here for uncorrelated components) of with constrained optimization <cit.>. 




§ RELATIONSHIP BETWEEN SCATTER MATRICES AND BETWEEN- AND WITHIN-SUBJECT COVARIANCE MATRICES



Comparing Eq. (<ref>) to (<ref>) yields:


    R_T = N  S_B .
 

One can also relate R_W to S_T as follows.  Expanding the definition (<ref>) yields (with abbreviated notation for the sums):


    R_W    =   ∑_il x_i^l  x_i^l ⊤ 
    - ∑_il x_i^l x̅_*^l ⊤  
    - ∑_ilx̅_*^l  x_i^l ⊤
    + ∑_ilx̅_*^l x̅_*^l ⊤
       =   ∑_il x_i^l  x_i^l ⊤ 
    - T ∑_ lx̅_*^l x̅_*^l ⊤
       =   ∑_il x_i^l  x_i^l ⊤ 
    - T ∑_ lx̅_*^l x̅_*^l ⊤  
    + T ∑_ lx̅_*^* x̅_*^* ⊤  
    - T ∑_ lx̅_*^* x̅_*^* ⊤
       =   ( 
        ∑_il x_i^l  x_i^l ⊤ 
    + T ∑_ lx̅_*^* x̅_*^* ⊤) 
    -
    ( 
        ∑_ilx̅_*^l x̅_*^l ⊤  
    - T ∑_ lx̅_*^* x̅_*^* ⊤)
    
       =   ∑_il ( x_i^l - x̅_*^*)
                  ( x_i^l - x̅_*^*)^⊤ 
    -
        ∑_il (x̅_*^l - x̅_*^*)
                  (x̅_*^l - x̅_*^*)^⊤
       =    S_T  - T  S_M .
 

Matrix S_M, defined in (<ref>), captures the variance of the mean across subjects and is zero when all subjects have equal mean: x̅_*^l = x̅_*^* .[ 
For time domain signals, this is often the case as it is customary to subtract the mean value (constant offset) to make the signals zero-mean. For rating data, an unequal mean indicates rater bias, i.e. a rater gives systematically lower/higher ratings than other raters, and it may be advantageous to remove this bias by equalizing all mean values. Note that, when S_M=0, the present definition of ISC is equivalent to the classic definition of intra-class correlation (ICC, see Appendix B). But in general, this matrix is not zero, which is important when generalizing the approach to non-linear transformations. There, enforcing zero mean in x does not guarantee zero mean in y.] 

With this result, scatter and covariance matrices can be related as follows:


    S_B    =   1/N R_T = 1/N ( R_W +  R_B) , 
     S_W    =    S_T -  S_B  
    =  R_W + T  S_M -  S_B 
    
          
    = N-1/N R_W - 1/N R _B   + T  S_M .






§ REGULARIZATION AND ROBUSTNESS



Note that a high inter-subject correlation corresponds to large eigenvalues of 𝐑_W^-1𝐑_B. These eigenvalues will be large if the selected direction exhibits large power in the space of 𝐑_B, or low power in the space of 𝐑_W.  The former is meaningful and desired, but the latter can occur with rank-deficient or impoverished data and may lead to spuriously high inter-subject correlation.   Truncating the eigenvalue spectrum of 𝐑_W has the desired effect of forcing the extracted data dimensions to have both high power in 𝐑_B and 𝐑_W.  This makes it less likely to find spuriously reliable dimensions.   Formally, the TSVD approach is to select 𝐕 according to:


    R_B  V =R̃_W  VΛ,


where the within-subject covariance is now regularized with the following approximation:


    R̃_W = Ũ_WΛ̃_WŨ_W^⊤,


where
Ũ_W=[  u_1  … u_K  ],  and 
Λ̃_W = 𝚍𝚒𝚊𝚐[ λ_1 …λ_K ]
are the K principal eigenvectors of R_W
and associated eigenvalues.    The inverse of R_W is computed according to:


    R̃_W ^-1= Ũ_WΛ̃_W^-1Ũ_W^⊤,


which then left-multiplies (<ref>) and yields the regularized solution to the CorrCA problem.  

The parameter K is the number of eigenvectors to retain in the representation of R_W.  Decreasing this number strengthens the level of regularization, and K=D is equivalent to not regularizing. 


Shrinkage regularization operates by a similar principle: dimensions of the data exhibiting low variance (corresponding to low eigenvalues of R_W) are enriched by incrementing these eigenvalues:  


    R̃_W =  U_WΛ̃_W U_W^⊤,

where
Λ̃_W= (1-γ) Λ_W + γλ̅ I ,
γ is the shrinkage parameter with 0 ≤γ≤ 1 , and λ̅ = Tr( R_W) / D is the mean eigenvalue of the (unregularized) within-subjects covariance. The effect of shrinkage, which retains the full dimensionality of the covariance matrix, is that the smaller eigenvalues are increased and larger eigenvalues are reduced (hence the term “shrinkage”).  

The advantage of shrinkage regularization is the simplicity of its implementation, where the regularized covariance can be simply computed as


    R̃_W = (1-γ)  R_W + γλ̅ I


without any further modifications to the generalized eigenvalue routine. TSVD instead requires more careful handling of the null-space when computing the general eigenvalues (see provided code). On the other hand, TSVD has the important advantage that the solutions V to the regularized eigenvalue problem (<ref>) diagonalize R̃_W as well as the original R_W. This means that the extracted components are strictly uncorrelated, whereas for shrinkage that is only approximately the case.

The performance of TSVD and shrinkage can be evaluated by maximizing ISC on one data set (training set) and using the resulting projection vectors to measure ISC on a different data set (test set). We do this on one of the application examples in section <ref>. 

Note that the TSVD approach can also be used to regularize the MCCA algorithm. In this case, the regularized inverse D̃ has to be computed separately for each block of the block-diagonal matrix D in (<ref>).





§ TESTING FOR STATISTICAL SIGNIFICANCE




The assessment of statistical significance of ρ requires knowledge of their distribution under the null hypothesis. The null hypothesis associated with CorrCA states that there exists no shared one-dimensional subspace in which the data of at least two subjects are correlated. Estimating the null distribution is complicated by two factors. First, CorrCA has free parameters, the optimization of which inflates the ISC values under the null hypothesis to an extent governed by the dimensionality of the problem (parameters T, D, and N). Second, temporal correlations in the data lead to larger variances than expected under the assumption of independent and identically distributed (IID) data typical for parametric statistical tests, and, thereby, to false positive detections of significant ISC values. The IID assumption may be reasonable for rating data, where rated individuals may be independent from one another, but it is not fulfilled for time series data such as EEG recordings. 
Here, we present three approaches to estimate the significance of individual correlated components as well as the dimensionality of the correlated subspace for IID as well as auto-correlated data.

Parametric null distribution—For IID data, we can use the F-statistics given in Eq. (<ref>) to test for significant ISC (provided the means are equalized). Evidently, this statistic is not applicable for ρ values obtained on training data, as the we have adapted parameters to optimize this statistic itself. However, we can use the optimal projection vectors to assess statistical significance on separate, previously unseen test data. An important limitation is that the T exemplars have to be independent. This is true perhaps for rating data, but typically not true for temporal signals, where samples in time are often correlated. For such data, the non-parametric test discussed above are more appropriate. To avoid overestimating the number of components due to multiple testing (there are always D components to be tested), we adopt a Bonferroni correction. This  limits the family-wise error rate, i.e. the probability of making one or more false discoveries among the D tests. Thus, only S values with p<α/D are considered significant, where α is the desired significance level and p is computed using the F-statistics.

Phase-scrambled surrogate data— For auto-correlated data we have to use shuffle statistics. This approach was introduced by <cit.> as a tool to test for nonlinearities in time series data.
The idea is to create data that are consistent with the null hypothesis but otherwise resemble the observed experimental data in terms of temporal and spatial correlations. The original surrogate data of <cit.> preserves the amplitude spectrum of the original data but uses a random phase spectrum, which is achieved through an application of the Fourier transform and its inverse. This approach was extended to additionally maintain the spatial covariance structure of multivariate data by using identical random phase shifts for all variables by <cit.>. Here, we use the approach of  to test for significant ρ based on the consideration that no inter-subject correlation can be present after phase-scrambling the data of each subject. Variants of this approach have previously been used to test for significant correlations <cit.>. The p-value for the ISC of a given correlated component is defined as the frequency with which the ISC computed from the random-phase surrogate data is exceeded by the maximum ISC of the most strongly correlated component in the
original data. Note that for every shuffle we obtain a set of ISC values, ρ_d, for the dth component. When we compare the ISC value on the original data, we measure p-value based on the distribution of maximum ρ, i.e. all ρ_d values are compared against the Null distribution of ρ_1. With this approach no further correction for multiple comparison is needed. 


Random circular shifts—A simpler approach to create surrogate data for auto-correlated data is to circularly shift the samples of each subject by a different offset, where the same offset is used for all variables of the same subject (i.e. sample index i=1 … T are shifted by a random offset o: i ← (i+ o) mod  T + 1). Just as the phase-scrambling approach, this procedure ensures that spatial correlations and spectral properties of the original data are maintained. Additionally, non-linear properties of the original data are maintained. This approach has previously been used in <cit.>. P-values are defined analogous to the phase-scrambling approach.

We test the validity of the three methods on simulated data in section <ref>. The goal is to determine under which conditions these methods determine the correct number of underlying correlated components.
 






§ VALIDATION ON SIMULATED DATA: DATA GENERATION AND EVALUATION METHODS



The generated signal and noise components, s^l(t) ∈ℝ^K and n^l(t) ∈ℝ^D, t = 1, , T, l=1, , N, were mapped to the measurement space as x_s^l(t) =  A_s^l  s^l(t), x_n^l(t) =  A_n^l  n^l(t) and normalized as x_s^l(t) ← x_s^l(t) / || x_s^l||_F, x_n^l(t) ← x_n^l(t) / || x_n^l||_F, where || x||_F denotes the Frobenius norm of the multivariate data set x(t). The mixing matrices A_s^l =  O_s^l  D_s^l ∈ℝ^D × K and A_n^l =  O_n^l  D_n^l ∈ℝ^D × D were generated as products of matrices O_s^l ∈ℝ^D × K (O_n^l ∈ℝ^D × D) with random orthonormal columns and diagonal matrices D_s^l ∈ℝ^K × K (D_n^l ∈ℝ^D × D), the diagonal entries of which were drawn randomly as D_ii = exp (d_i), d_i ∼𝒩(0, 1) and normalized to a maximum of 1. This ensured that the non-zero eigenvalues of the signal and noise covariances matrices decayed exponentially as is the case for many real-world data sets. Unless otherwise noted, the same mixing matrices A_s^l ≡ A_s and A_n^l ≡ A_n were used for all subjects, reflecting the fundamental assumption of CorrCA that the spatial covariance structure of signal and noise components is the same in all subjects. The emulated data were generated as x^l(t) = ξ x_s^l(t) + (1-ξ)  x_n^l(t), where the scalar parameter ξ = 10^SNR/20 / (1 + 10^SNR/20) was used to adjust the SNR. Note that SNR was defined here in terms of the contribution of the signal and noise components to the Frobenius norm of the entire multivariate measurement. It is, therefore, different from the definition provided in Eqs. (<ref>) and (<ref>), which applies to single channels or components and relates to the ISC through Eq. (<ref>). Note that in most of our experiments, the ISC of each simulated signal component was set to 1, corresponding to infinite SNR according to Eq. (<ref>).

We applied CorrCA to the emulated data and measured its performance in terms of the achieved ISC (average over the first K components) on the training set. We also measured the average ISC for a test set, i.e., using the same CorrCA projections V on new emulated data of identical size. Moreover, we measured the performance of CorrCA to reconstruct the true (emulated) correlated components s^l(t) and their corresponding mixing matrices A_s^l by comparing them to y^l(t) and the estimated forward model A (section <ref>). The performance metric for this was the Pearson correlation between individual simulated and reconstructed components or columns of A_s. With this, we measure if each component was correctly identified. However, because components with equal ISC can be arbitrarily mixed within a subspace (Appendix <ref>), we also assessed identifiability in terms of the angle between the subspaces spanned by the first K columns of the estimated Â and A_n, as well as the angle between the simulated components s^l(t) and the first K dimensions of y^l(t) reconstructed by CorrCA. Angles were normalized to the interval [0, 1]. Finally, we estimated the number of simulated correlated components, K, using phase-scrambled surrogate data, surrogate data obtained using random circular shifts, and the parametric F-test. Empirical null distributions were obtained using 1000 surrogate data sets. The number of components with p-values smaller than α = 0.05 was used as an estimate for K. For the parametric test, the T samples were randomly split into training and test sets of equal size, where the CorrCA projections V were obtained on the training set, and the statistical assessment was conducted on the test set. For this approach, the median number of components with p-values smaller than α = 0.05 across 100 random train/test splits was used as an estimate of K.

All experiments were repeated 100 times. Mean values and standard deviations of all performance measures are depicted in Figures <ref>—<ref>.




§ EXTENSION TO NON-LINEAR MAPS USING KERNELS



So far we have only considered linear transformations between x and y. Now, given the close relationship between CorrCA and LDA, we can extend CorrCA to non-linear transformations following the approach of kernel-LDA <cit.>. Assume a non-linear mapping Φ(·) for the data of subject l:


    Z^l = Φ( X^l) .


In this notation, X^l has dimension D × T spanning the entire multivariate time series for subject l, and Z^l has dimensions D' × T'. Thus, the non-linear transformation Φ(·) is a mapping from space-time to a new space with possibly quite different dimensions. 
[In the example of brain signals, Z^l could represent, for instance, the power spectrum in source space <cit.>, where D' captures the number of sources and T' the number of frequency bins. In the application examples and in the code we simplify this by applying the non-linear mapping to each time point separately, so that T=T' and we have a purely spatial nonlinear mapping, similar to conventional kernel-LDA. However, this more general formulation allows one in principle to reduce the time dimension as well. This is important because the computational complexity of the kernel approach scales with O(T^3), and thus becomes prohibitive for long time-domain signals.] 
The approach now is to first extract non-linear features of the data with Φ(·), and then combine these features linearly with projection v: 


    y^l =  v^⊤ Z^l.


Note that the transformation Φ(·) does not need to be defined explicitly.  Instead, the “kernel trick” <cit.> allows one to specify this non-linearity only implicitly, by defining instead the inner product of vectors in this new D'-dimensional feature space with a kernel function K( · , · ):


    Z^k^⊤ Z^l =  K( X^k,  X^l) .
  

Knowing how to compute the inner product is sufficient, if all the expressions in the algorithm can be formulated in terms of this inner product.[Note that this inner product is only over dimension D'. The function K( · , · ) is therefore only a "kernel" in this D' dimensional space.] In a slight abuse of notation, we will write K^kl =  K( X^k,  X^l) as shorthand for this matrix defined for each pair of subjects lk. Note that the dimensions of K^kl are T' × T' and thus this matrix can potentially be quite large. In cases where the non-linear transformation Φ(·) can be expressed explicitly, and D' ≪ T', one may chose to operate directly with Z^l. Otherwise, it will be beneficial to leverage the definition of the non-linear transformation in terms of K^kl. In the following, we will rewrite the optimization criterion ρ in terms of K^kl alone.  

The crucial step is to define the projection vectors in terms of the samples of the mean across subjects. 

    v = ∑_i=1^T'α_i z̅^*_i = Z̅^* ,

where z̅^*_i is the across-subject mean vector for exemplar i in the non-linear mapped space (ith column of Z̅^*).  Here, α_i are parameters that indicate how exemplars z̅_i are to be combined to represent the projection vector v. The approach will be to find optimal vector  instead of vector v, and thus we need to rewrite the optimality criterion in term of .[As an alternative to (<ref>), we could have also defined the projection vector in terms of the samples of all subjects:


    v = ∑_l=1^N ∑_i=1^T'α_i^l  z^l_i  .


The derivations can all be repeated analogously and the resulting algorithm is the same, but with NT' degrees for freedom in . In the provided code, this is referred to as the “full model”, whereas equations derived here for the expansion in terms of the mean (<ref>) are referred to as the “mean model”. In the example provided in Figure <ref>, the numerical results are very similar for both models. The question as to whether the mean model is sufficient is empirical and comes down to whether there is enough diversity in the mean response to capture the individual variations across subjects. Mathematically, this implies that the expansion in terms of the mean is complete, i.e., the mean response matrix Z̅^* is full rank.
] 

Combining (<ref>) with (<ref>) one can now express the non-linear components y as a linear combination of K^lk as follows:


    y^l    =    v^⊤ Z^l 
    =   
    ^⊤1/N∑_k=1^N  Z^k^⊤ Z^l
    
       =   ^⊤1/N∑_k=1^N  K^kl
    =
    ^⊤K̅^*l ,


where the bar in K̅^*l, indicates that we are taking the mean, and the asterisk specifies over which index this mean is taken. It analogously follows that,

    v^⊤Z̅^* = ^⊤K̅^** .

Denote the columns of matrix Z^l and K̅^*l as z_i^l and k̅_i^*l. Then we can also write with this notation for the average:


    v^⊤z̅_i^* = ^⊤k̅_i^** .


The within-subject and total variance of the non-linear features z_i^l can thus be written as:


    r_W    =    v^⊤ R_W  v
    =
    ∑_l=1^N∑_i=1^T' v^⊤ 
    (  z_i^l - z̅_*^l )
    (  z_i^l - z̅_*^l )^⊤ v
       =   ^⊤ C_W
     , 
    
    r_T    =    v^⊤ R_T  v
    =
    N^2 ∑_i=1^T' v^⊤ 
    ( z̅_i^* - z̅_*^* )
    ( z̅_i^* - z̅_*^* )^⊤ v
       =   ^⊤ C_T
     ,


where we defined the within-subject and total covariance of k_i^*l as:


    C_W    =   ∑_l=1^N∑_i=1^T' 
    ( k̅_i^*l - k̅_*^*l )
    ( k̅_i^*l - k̅_*^*l )^⊤ , 
     C_T    =    
    N^2 ∑_i=1^T' 
    ( k̅_i^** - k̅_*^** )
    ( k̅_i^** - k̅_*^** )^⊤ .


Because of the symmetry of these expressions with the definitions of R_W and R_T, we have again, C_B =  C_T -  C_W, and can therefore write the inter-subject correlation in the non-linear space as a function of the new parameters  as follows:


    ρ = 
     v^⊤ R_B  v/ v^⊤ R_W  v 
    =
    ^⊤ C_B /^⊤ C_W .


This can be solved again as an eigenvalue problem but now with the within- and between-subject covariance matrices of the kernel vectors k̅_i^*l and the corresponding projection vector .  Note that a similar result is obtained for the kernel version of CCA, namely, the conventional CCA algorithm is applied to kernel matrices instead of the original data <cit.>.  



To demonstrate the method, we generated a simple 2-dimensional example where the shared dimension is non-linear, and this non-linear relationship has to be discovered from the data alone. As with the example for linear CorrCA (Figure <ref>), we have simulated N=2 “subjects”, each providing two signals (D=2), now with T=40 samples in time (Figure <ref> A & B). In order to test whether the algorithm finds the correct non-linear mapping we generate data with a known non-linear relationship. In this case the two subjects share the amplitude (distance from the origin in the 2-dimensional input space of x; see Figure <ref> C). However, the phase angle in this 2-dimensional plane is selected randomly at each time point and for each subject. We apply kernel-CorrCA using Gaussian kernels <cit.> keeping the time axis unchanged (T=T'; see also our companying code). The first component dimension found by kernel-CorrCA, y_1, is approximately linear with the amplitude and is independent of phase in the original 2D plane (Figure <ref> G & H). The algorithm therefore discovered in the first component amplitude as the shared non-linear dimension. The second component dimension, y_2 is essentially a nonlinear function of component y_1 (Figure <ref> F), hence it also has high ISC. The specific non-linearity between the two components is arbitrary, but does enforce that the two variables are uncorrelated. We obtain nearly identical results regardless of whether we use Gaussian or Tanh kernels <cit.>, although the specific relationship between y_1 and y_2 differs.  







§ FURTHER REMARKS ON INTER-SUBJECT CORRELATION




 §.§ Relationship of ISC to Intra-Class Correlation



The definition of inter-subject correlation in the present work is similar to the definition of intra-class correlation (ICC) first introduced by <cit.>. Fisher's ICC is referred to as pairwise inter-class correlation and differs from Pearson's correlation coefficient <cit.>. A more modern definition of intra-class correlation, which in the present context would allow for a different number of subjects per sample (i.e. N_i), was introduced by <cit.>:


    r_p = ∑_i=1^T1/N_i-1∑_k=1^N_i∑_l=1, k ≠ l^N_i (y_i^l - y̅_*^*)(y_i^k - y̅_*^*)/∑_i=1^T∑_l=1^N_i (y_i^l - y̅_*^*)^2.


The definition of pairwise ICC, r_p, differs from our definition of ISC, ρ, in that it subtracts the total  mean y̅_*^*, whereas we subtract the sample mean of each subject, y̅_*^l. The two are identical in the case of unbiased ratings, i.e. S_M=0 defined in Eq. (<ref>). A definition comparable to ISC can be found in <cit.> in the context of MCCA under the term “modified inter-class correlation”, but is not commonly seen in the literature. Instead, the generalized eigenvalue solution to MCCA presented here has been often derived using the sum of correlations, subject to a normalization constraint <cit.>.

In Table <ref>  we summarize various uses of these correlation measures. For each case the indices play different roles. In the context of ICC, index i represents different classes whereas l,k are different exemplars of a class. Each class can have a different number of exemplars N_i. In the context of ISC, l, k are subjects and i represent measures taken, perhaps in time. In the case of ratings, l,k are the raters and i are the items being rated. In <cit.>, the measure r_p is used to assess intra-familial correlations (correlations between siblings). In this case, k,l indexes the siblings, whereas i enumerates the families. The algorithms listed in Table <ref> aim to extract components in multi-dimensional data that maximize the corresponding statistic. Where no algorithm is listed, either CorrCA or LDA could be used. 

 




The definition of ISC/IRC in Eq. (<ref>) can be thought of as a Pearson correlation coefficient, where we concatenate the signals (after subtracting the individual subject/rater means) corresponding to all possible pairs of subjects. In contrast to conventional Pearson correlation, with this symmetric definition (where we include pairs lk as well as kl), the signals of different subjects have to have identical scale in order to achieve perfect ISC/IRC (see section <ref>). In the case of ratings, this means that different raters can only differ in their mean ratings but have to otherwise provide identical deviations from the mean in order to achieve perfect inter-rater correlation. This means that this metric is sensitive to multiplicative bias (i.e. each rater/subject has a different scale). The issue of bias in the context of inter-rater agreement is discussed in more detail by <cit.>, who argue that Pearson correlation should be used to assess test-retest reliability to account for learning effects of subjects from one test to the next. If there is an undesirable proportional bias, then one can easily correct for this bias by standardizing the data, i.e. dividing by the standard deviation for each subject/rater in each dimension prior to applying CorrCA. 

As we have shown here, maximizing IRC is equivalent to maximizing inter-rater-reliability, where we have defined reliability in Eq. (<ref>) as the variance of the mean over the mean of the variance. For instance, if the rating is the speed of performing a task (as in the application <ref>), this would mean that a 100 millisecond variation across repeats is fairly reliable for a task that takes 5 seconds to complete, but really quite unreliable if the task only takes 0.5 second in average. 

It should be noted that, when ratings are on a discrete numerical scale, IRC may not be an appropriate measure of reliability. For instance, a rating on an integer scale of 1, 2, 3, 4, 5, may have an inter-rater variability of ± 1.0. With the present definition of IRC — ρ in Eq. (<ref>) — this same variability would be considered half as reliable if the mean rating is 4 as compared to a mean of 2. Whether this is desirable or not may depend on the specific ratings statistic. Finally, IRC can only be used meaningfully for ratings that are on a proportional scale. For categorical ratings one should use instead the conventional Cohen's κ or similar inter-reliability measures.



 §.§ ISC is normalized



To our knowledge, the present definition of ρ in Eq. (<ref>) does not appear in this exact form in the literature, except for our previous work <cit.>. Therefore, we take some time here to show that this definition of correlation has indeed a maximal value of 1, and that this occurs only when the all signals are identical across subjects (not just proportional as in Pearson's correlations), except for a subject-dependent mean, which is permitted.

Without loss of generality, let us assume that all signals are zero mean, y̅^l_*=0. To show that ρ≤ 1, we can also demonstrate, equivalently, that the following expression is non-negative:


    0    ≤     (N-1) r_W - r_B  
    
       =    
    (N-1) ∑_l ∑_i  (y^l_i)^2 - ∑_l ∑_k≠ l∑_i y^l_i y^k_i 
    
       =    
    N ∑_l ∑_i  (y^l_i)^2 - ∑_i  ( ∑_l  y^l_i )^2 =E .


The last expression is abbreviated as E. To find the y^l_i that minimizes E we find the location where the gradient is zero and check whether the curvature at that location is non-negative. 


    ∂ E/∂ y_j^k   =    
    2 N y^k_j - 2 ∑_i   ∑_l  y^l_i ∑_l δ^lk_ij
       =    
    2 N y^k_j - 2 ∑_l  y^l_i   ,
    
    ∂^2 E/∂ y_j^k ∂ y_j^l   =    
    2Nδ^kl_ji - 2 ∑_n δ^kn_ij
       =    
     2δ_ij (Nδ^kl-1).


Solving for ∂ E/∂ y=0  yields an unique solution:


    y_j^k = 1/N∑_l y_j^l = y̅^*_j.


This solution is a linear manifold, as one can pick arbitrary scale for y̅^*_j. Inserting this into (<ref>) yields E=0 at that point. To determine whether this is indeed a minimum value for E, and hence ρ≤ 1, we need to show that the Jacobian (<ref>) has only non-negative eigenvalues. Because of the δ_ij in (<ref>), we can find the eigenvalues for each ij separately, i.e. the eigenvalues of J = N   I -  1, where I is the N-dimensional identity matrix and 1 is a N-dimensional vector with all values set to 1. This J is the same for all ij, and the eigenvalue equation is for y∈ℝ^N: 


    J y = λ y .


The characteristic polynomial for these eigenvalues reads:


    0   =   det( N I -  1 1^⊤ - λ I)) 
    
       =   
    (N-λ)^N det(  I - (N-λ)^-1 1 1^⊤) 
    
       =   
    (N-λ)^N ( 1 - (N-λ)^-1 N) 
    
       =   - (N-λ)^N-1λ .
 

This means that there is a (N-1)-fold eigenvalue, λ=N, and a unique eigenvalue, λ=0. Therefore, the Jacobian is semidefinite, i.e. it has positive curvature in all directions, except in the direction of the eigenvectors of λ=0. Inserting this in (<ref>) yields y=c  1 with a free parameter c. If we compare that to (<ref>), we see that these are the same, with c=y_i^*. This means that, when all signals are the same for all subjects, then we achieve a global minimum for E and, therefore, ρ is maximal and equals to 1. The curvature is zero in one direction at the minimum because one can change the global scale of all signals without changing E=0.   



 §.§ ISC between individuals and a group of subjects


 
In a number of studies, we compare brain responses of individual subjects to that of the group by correlating between pairs of subjects, and then averaging over pairs involving the subjects of interest. This correlation of individual subject with a group is then used to compare with individual measures <cit.>. For this purpose, we used the following definition for single-subject ISC <cit.>: 


    ρ_k    =   ∑_l=1,l  k^N (r_kl+r_lk)/∑_l=1,l  k^N (r_ll+r_kk) ,
    
    r_kl   =   ∑_i=1^T (y_i^k - y̅_*^k)(y_i^l - y̅_*^l).


This can be computed for each extracted component y_j(t)= v^⊤_j  x(t).



