Title: MEASURING CAUSAL INFLUENCE WITH BACK-TO-BACK REGRESSION: THE LINEAR CASE

Abstract: Identifying causes from observations can be particularly challenging when i) potential factors are difficult to manipulate individually and ii) observations are complex and multi-dimensional. To address this issue, we introduce "Back-to-Back" regression (B2B), a method designed to efficiently measure, from a set of co-varying factors, the causal influences that most plausibly account for multidimensional observations. After proving the consistency of B2B and its links to other linear approaches, we show that our method outperforms least-squares regression and cross-decomposition techniques (e.g. canonical correlation analysis and partial least squares) on causal identification. Finally, we apply B2B to neuroimaging recordings of 102 subjects reading word sequences. The results show that the early and late brain representations, caused by low-and high-level word features respectively, are more reliably detected with B2B than with other standard techniques.

1. INTRODUCTION:

Natural sciences are tasked to find, from a set of hypothetical factors, the minimal subset that suffices to reliably predict novel observations. This endeavor is impeded by two major challenges.

First, causal and non-causal factors may be numerous and collinear. This issue becomes increasingly pronounced as the number of potential factors increases. In neuroscience, for example, identifying whether the frequency of a word presented on a subject's retina modulates brain activity can be surprisingly difficult. Indeed, the frequency of words in natural language covaries with other factors such as their length (short words are more frequent than long words), their categories (determinants are more frequent than adverbs) and so forth (Kutas & Federmeier, 2011; Pegado et al., 2014) . Instead of selecting a set of words that control for all of these factors simultaneously, it is thus common to use forward modeling, i.e. to train a model to predict observations (e.g. brain activity) from a minimal combination of competing factors (e.g. word length, word frequency, e.g. (Huth et al., 2016) ), and investigate, in the model, the estimated contribution of each factor (Friston et al., 1994) .

The second challenge to measuring causal influence is that observations can be large and complex. The relationship between causes and effects is thus often considered in a backward manner, by training models to maximally predict causes from multidimensional observations. For example, brain activity is often recorded with hundreds or thousands of sensors simultaneously. As multiple sensors may be affected by common noise sources, it is common to use backward regression, by, for example, fitting a support vector machine across multiple sensors to decode the category of a stimulus (Cichy et al., 2014; Kriegeskorte et al., 2008; Norman et al., 2006) .

Both forward and backward modeling have competing benefits and drawbacks. Specifically, forward modeling disentangles the independent contribution of collinear factors, but does not combine multidimensional observations. By contrast, backward modeling combines multiple observations, but does not disentangle collinear factors (Weichwald et al., 2015; Hebart & Baker, 2018; King et al., 2018) . To combine the benefits of forward and backward modeling, several authors have proposed to use cross-decomposition techniques such as Partial Least Squares (PLS) and Canonical Correlation Analysis (CCA) (de Cheveigne et al., 2019) . CCA and PLS aim to find, from two sets of data X and Y , the components H and G were XH and Y G are maximially correlated or maximally covarying respectively. Because CCA and PLS are based on a generalized eigen decomposition, their resulting coefficients are mixing the features of X and Y in a way that makes them notoriously difficult to interpret (Lebart et al., 1995) .

= × + × Y ∈ R n×dy X ∈ R n×dx E ∈ D dx×dx N ∈ R n×dx F ∈ R dx×dy Observations Factors Cause selection Noise Cause-effect mapping X ∼ P (X) N ∼ P (N ) Ê = diag((X 2 X 2 + Λ X ) -1 X 2 Y 2 1) X: regression from Y to X (Y 1 Y 1 + Λ Y ) -1 Y 1 X 1 2) regression from X to X )

Figure 1 : Back-to-back regression identifies the subset of factors E ii = 1 in X that influence some observations Y by 1) regressing from Y to X to obtain X, and 2) returning the diagonal of the regression coefficients from X to X.

Here, we introduce the 'back-to-back regression' (B2B), which not only combines the benefits of forward and backward modeling (Section 2), but also provide robust, interpretable and unidimensional causal coefficients for each tested factor.

The present paper focuses on the restricted issue of disentangling the causal influence of linearly correlated predictors (X) onto noisy multivariate observations (Y ). The present approach thus differs from other causal discovery algorithms based on temporal-delays and/or nonlinear interactions in systems where the directionality of causation (from X to Y or vice versa) is unknown (e.g. (Peters et al., 2017; Granger, 1969; Janzing et al., 2013; Schölkopf et al., 2016) .

After detailing B2B method and proving its convergence (Section 2.2), we show with synthetic data that it outperforms state-of-the-art forward, backward and cross-decomposition techniques in identifying causal influence (Section 3.1). Finally, we apply B2B to a large neuroimaging dataset and reveal that distinct but collinear word features lead to distinguishable brain representations (Section 3.5).

2. BACK-TO-BACK REGRESSION:

We consider the measurement of multivariate signal Y ∈ R n×dy , generated from a set of putative causes X ∈ R n×dx , via some unknown linear apparatus F ∈ R dx×dy . Not all the variables in X exert a causal influence on Y . By considering a square binary diagonal matrix of causal influences E ∈ D dx×dx , we denote by XE the causal factors of Y . In summary, the problem can be formalized as:

EQUATION

)

where i is a given sample, and n i is a sample-specific noise drawen from a centered distribution. While the triplet of variables X and N are independent, we allow each of them to have any form of covariance. In practice, we observe n samples (X, Y ) from the model. This problem space, along with the sizes of all variables involved, is illustrated in Figure 1 . Given the model in Equation (1), the goal of Back-to-Back Regression (B2B) is to estimate the matrix of causal influences E.

2.1. ALGORITHM:

Back-to-Back Regression (B2B) consists of two steps. First, we estimate the linear regression coefficients Ĝ from Y to X, and construct the predictions X = Y Ĝ. This backward regression recovers the correlations between Y and each factor of X. Second, we estimate the linear regression coefficients Ĥ from X to X. The diagonal of the regression coefficients Ĥ, denoted by Ê = diag( Ĥ), is the desired estimate of the causal influence matrix E, as detailed in the A.1.

If using l2-regularized least-squares (Hoerl, 1959; Rifkin & Lippert, 2007) , B2B has a closed form solution:

EQUATION

where Λ X and Λ Y are two diagonal matrices of regularization parameters, useful to invert the covariance matrices of X and Y if these are ill-conditioned.

Performing two regressions over the same data sample can result in overfitting, as spurious correlations in the data absorbed by the first regression will be leveraged by the second one. To avoid this issue, we split our sample (X, Y ) into two splits (X 1 , Y 1 ) and (X 2 , Y 2 ). Then, the first regression is performed using (X 1 , Y 1 ), and the second regression is performed using (X 2 , Y 2 ). To compensate for the reduction in sample size caused by the split, B2B is repeated over many random splits, and the final estimate Ê of the causal influence matrix is the average over the estimates associated to each split (Breiman, 1996) . To accelerate this ensembling procedure, we implemented an efficient leave-one-out cross-validation scheme as detailed in (Rifkin & Lippert, 2007) , as follows:

Ŷ = (Σ X GY -diag(Σ X G)Y ) / diag(I -Σ X G) (element-wise division) (4)

where Σ X is the X kernel matrix and where G is computed with an eigen decomposition of X:

Σ X = QV Q T G = Q(V + λI) -1 Q T (5)

where Q, V and λ are the eigen vectors, eigen values and regularization, respectively.

We summarize the B2B procedure in Algorithm 1. The rest of this section provides a theoretical guarantee on the correctness of B2B.

Algorithm 1: Back-to-back regression.

Input: input data X ∈ R n×dx , output data Y ∈ R n×dy , number of repetitions m ∈ N. Output: estimate of causal influences Ê ∈ D dx×dx . 1 Ê ← 0; 2 for i = 1, . . . , m do 3 (X, Y ) ← ShuffleRows((X, Y )); 4 (X 1 , Y 1 ), (X 2 , Y 2 ) ← SplitRowsInHalf((X, Y )); 5 Ĝ = LinearRegression(Y 1 , X 1 ) ; Ĝ = (Y 1 Y 1 + Λ Y ) -1 Y 1 X 1 6 Ĥ = LinearRegression(X 2 , Y 2 Ĝ) ; Ĥ = (X 2 X 2 + Λ X ) -1 X 2 Y 2 Ĝ 7 Ê ← Ê + diag( Ĥ); 8 end 9 Ê ← Ê/m; 10 Ŵ ← LinearRegression(X Ê, Y ); 11 return Ê, Ŵ

2.2. THEORETICAL GUARANTEES:

Theorem 1 (B2B consistency -general case). Consider the B2B model from Equation Y = (XE + N )F , N centred and full rank noise. Let Img(M ) refers to the image of the matrix M . If F and X are full-rank on the Img(E), then, the solution of B2B, Ĥ, will minimize min

EQUATION

Assuming, without loss of generality, that the active features in E are the k ∈ Z : k ∈ [0, d x ] first features, and rewriting X = (X 1 , X 2 ) and N = (N 1 , N 2 ) (X 1 and N 1 containing the k first features), we have:

EQUATION

where Σ AB is the covariance of A and B, and:

EQUATION

In the absence of noise, we have Σ N1N1 = 0, and so diag k ( Ĥ) = I, and

diag( Ĥ) = diag(E)

Therefore, we recover E from Ĥ.

In the presence of noise, the causal factors of E correspond to the positive elements of diag( Ĥ). The methods to recover them are presented in the Appendix A.4.

3. EXPERIMENTS:

We perform two sets of experiments to evaluate B2B: one on controlled synthetic data, and a second one on a real, large-scale magneto-encephalography (MEG) dataset. We use scikit-learn's PLS and RidgeCV (Pedregosa et al., 2011) as well as pyrcca's regularized canonical component analysis (RegCCA, (Bilenko & Gallant, 2016 )) objects to compare B2B against the standard baselines.

3.1. SYNTHETIC DATA:

We evaluate the performance of B2B throughout a series of experiments on controlled synthetic data. The purpose of these experiments is to evaluate the ability of B2B in terms of prediction of independent and identically distributed data, as well as a method to recover causal factors.

The data generating process for each experiment constructs n = 1000 training examples according to the model Y = (hXE + N )F , where h is a scalar that modulates the signal-to-noise ratio.

Here, F ∈ R dx×dy contains entries drawn from N (0, σ 2 ) where σ 2 is inversely proportional to

d x , X ∈ R n×dx contains rows drawn from N (0, Σ X ), N ∈ R n×dx contains rows drawn from N (0, Σ N ), E ∈ R dx×dx is a binary diagonal matrix containing n c ones, Σ X = AA where A ∈ R dx×dx contains entries drawn from N (0, σ 2 ), Σ N = BB where B ∈ R dx×dx contains entries drawn from N (0, σ 2 )

, and the factor h ∈ R + .

To simulate a wide range of experimental conditions, we sample 10 values in log-space for

d x , d y ∈ [10, 100], n c ∈ [3, 63], h ∈ [0.001, 10].

We discard the cases where n c > d x , limit d x , d y to 100 to keep the running time under 2 hours for each condition, and average over 5 random seeds.

We compare the performance of B2B against four competing methods, all implemented in scikit-learn (Pedregosa et al., 2011) and pyrcca (Bilenko & Gallant, 2016):

3.2. BASELINE MODELS:

Forward regression consists of an l2-regularized "ridge" regression from the putative causes X to the observations Y : Backward regression consists of an l2-regularized "ridge" regression from Y to X:

H f wd = (X T X + λI) -1 X T Y (10)

EQUATION

CCA finds G cca ∈ R dz,dy and H cca ∈ R dz,dx s.t. X and Y are maximally correlated in a latent Z space:

EQUATION

PLS finds G pls ∈ R dz,dy and H pls ∈ R dz,dx s.t. X and Y are maximally covarying in a latent Z space:

EQUATION

We employ 5-fold cross-validation to select the optimal number of components for CCA and PLS.

Regressions were 2-regularized with a λ regularization parameters fitted with the efficient leave-oneout procedure implemented in scikit-learn RidgeCV (Pedregosa et al., 2011) .

3.3. EVALUATING CAUSAL DISCOVERY FROM MODELS' COEFFICIENTS:

B2B leads to unbiased (i.e. zeros-centered) scalar coefficients for non-causal features. In contrast, the Forward, Backward, CCA and PLS models lead to a loading vector H i per feature i (or one vector G i for the backward model). To transform such vector into an estimated causal contribution Ê, we take the sum of square coefficients:

Êi = j H j i 2

To estimate whether models accurately identify causal factors, we compute the area-under-the-curve (AUC) across factors AU C(E, Ê). The AUC allows evaluating the capacity of models at detecting the causal importance of factors when ground truth labels are available, as is the case in this setup.

We report AUC results in Figures 2 (top) and 5 (left, in Appendix), and compare favorably to all baselines.

3.4. EVALUATING CAUSAL DISCOVERY WITH HELD-OUT PREDICTION RELIABILITY:

In most cases, E is not known and AUC can thus not be estimated. To address this issue, we assess the ability of each model to reliably predict independent and identically distributed data from Y , given all of the X features versus all-but-ones feature X -i (i.e. 'knock-out X'). This procedure Top. Average brain response to words (word onset at t=0 ms), as viewed from above the head (red= higher gradient of magnetic flux). Bottom. Each line represents a magnetometer, color-coded by its spatial position. Posterior responses, typical of primary visual cortex activity, peak around 100 ms after word onset and are followed by an anterior propagation of activity typical of semantic processing in the associative cortices.

results in two correlation metrics R f ull and R knockout , whose difference ∆R i = R f ull -R knockout indicates how much each X i improves the prediction of Y . In our figures, ∆R is the average of ∆R i . A higher score means that for prediction, the model relies on individual features rather than combinations of features.

We show in Appendix A.3 pseudo-code to assess feature importance for our algorithm as well as baselines. For the Backward Model, feature importance cannot be assessed as the X colinearity is never taken into account.

We show in Figures 2 (bottom) and 5 (right, in Appendix) that our method outperforms baselines.

3.5. MAGNETOENCEPHALOGRAPHY DATA:

Next, we apply our method to brain imaging data from the anonymized multimodal neuroimaging "Mother Of all Unification Studies" (MOUS) dataset (Schoffelen et al., 2019) . The dataset contains magneto-encephalography (MEG) recordings of 102 healthy native-Dutch adults who participated in a reading task. Twelve subjects were excluded from the analysis because of corrupted file headers. Subjects were exposed to a rapid serial visual presentation of Dutch words. The word lists consisted of 120 sentences, and scrambled lists of the same words. Each word was presented on the computer screen for 351ms on average (min: 300ms, max: 1400ms). Successive words were separated by a blank screen for 300ms, and successive sentences were separated by an empty screen for a few (3-4) seconds.

3.5.1. MEG PREPROCESSING:

The raw MEG data was bandpass-filtered between 0.1 and 40Hz using MNE-Python default parameters (Gramfort et al., 2013; 2014) . Specifically, we used a zero-phase finite impulse response filter (FIR) with a Hamming window and with transition bands of 0.1Hz and 10Hz for the low and high cut-off frequencies. The raw data was then segmented 100ms before word onset and 1s after word onset (t = 0ms corresponds to word onset). Finally, each resulting segment was baseline-corrected between -100ms and 0ms, and decimated by 5 and thus led a sampling frequency of 240Hz. The average responses across words is displayed in Figure 3 . For each subject and each time sample relative to word onset, we build an observation matrix Y ∈ R n×dy of n ≈ 2700 words by d y = 301 MEG channels (273 magnetometers and 28 compensation channels). Each of the columns of Y is normalized to have zero mean and unit variance.

3.5.2. FEATURE DEFINITION:

We aim to identify the word features that cause a variation in brain responses. We consider four distinct but colinear features. First, 'Word Length' refers to the total number of letters. Word Length is expected to specifically cause a variation in the early evoked MEG responses (i.e. from 100 ms after stimulus onset) elicited by the retinotopically-tuned visual cortices (e.g. (Pegado et al., 2014) .). Second, 'Word Frequency' indexes how frequently each word appears in Dutch and was derived with the the Zipf logarithmic scale of (Van Heuven et al., 2014) provided by the WordFreq package (Speer et al., 2018) . Word Frequency is expected to specifically cause a variation in the late evoked MEG responses (i.e. from 400 ms), because it variably engages semantic processes in the temporal cortices (Kutas & Federmeier, 2011) . Third, 'Word Function' indicates whether each word is a content word (i.e. a noun, a verb, an adjective or an adverb) or a function word (i.e. a preposition, a conjunction, a determinant, a pronoun or a numeral), and was derived from Spacy's part of speech tagger (Honnibal & Montani, 2017) . To our knowledge, this feature has not been thouroughly investigated with MEG. Its causal contribution to reading processes in the brain thus remains unclear. Finally, to verify that B2B and other methods would not inadequately identify non-causal features, we added a dummy feature, constructed from a noisy combination of Word Length and Word Frequency: dummy = z(length) + z(f requency) + N , where z normalizes features and N is a random vector sampling Gaussian distribution (all terms thus have a zero-mean and a unit-variance). This procedure yields an X ∈ R n×dx matrix of n ≈ 2700 words by d x = 4 features for each subject. Each of the columns of X is normalized to have a mean and a standard deviation of 0 and 1 respectively.

3.5.3. MODELS AND STATISTICS:

We compare B2B to four standard methods: Forward regression, Backward regression, CCA and PLS, as implemented in scikit-learn (Pedregosa et al., 2011) and (Bilenko & Gallant, 2016) , and optimized with nested cross-validation over twenty l2 regularization parameters logarithmically spaced between 10 -4 and 10 4 (for regression and CCA methods) or 1 to 4 canonical components (for PLS).

We used the feature importance described in Algorithm 2 to assess the extent to which each feature X i specifically improves the prediction of held-out Y data, using a 5-fold cross-validation (with shuffled trials to homogeneize the distributions between the training and testing splits).

Each model was implemented for each subject and each time sample independently. Pairwise comparison between models were performed using a two-sided Wilcoxon test across subjects (n=90) using the average ∆R across time. Corresponding effect sizes are shown in Figure 4 , and p-values are reported below.

3.5.4. RESULTS:

We compared the ability of Forward regression, Backward regression, CCA, PLS and B2B to estimate the causal contribution of four distinct but collinear features on brain evoked responses to words.

As expected, the Backward model reveals a similar decoding time course for Word Length and Word Frequency, even though these features are known to specifically influence early and late MEG responses respectively (Kutas & Federmeier, 2011) . In addition, the same decoding time course was observed for the dummy variable. These results illustrate that backward modelling cannot be used to estimate the causal contribution of collinear features.

We thus focus on the four remaining methods (i.e. Forward Regression, PLS, CCA, and B2B) and estimate their ∆R (i.e. the improvement of Y prediction induced by the introduction of a given feature into the model,as described in Algorithm 2). Contrary to the Backward Model, none of the models predicted the Dummy Variable to improve the Y prediction: all ∆R < 0 (all p > .089).

Figure 4 shows, for each model, the effects obtained across time (left) and subjects (right).

Word Length and Word Frequency improved the prediction performance of all methods: ∆R > 0 for all models (all p < 0.0001). As expected, the time course associated with Word Length and Word Frequency rose from ≈ 100 ms and from ≈ 400 ms respectively. Furthermore, Word Function improved the prediction performance of all models (all p < 0.0002) except for PLS (p = 0.7989). Overall, these results confirm that Word Length, Word Frequency and Word Function causally influence specific periods of brain responses to words.

To assess which model would be most sensitive to these causal discoveries, we compared B2B to other models across subjects (Figure 4 right). For Word Length B2B outperforms all models (all p < 0.00001) but CCA (p = 0.0678). For Word Frequency, B2B outperforms all models (all p < 0.0006). For "Word Function", B2B outperforms all models (all p < 0.0015). Overall, these results show that B2B reliably outperforms standard methods, especially when the effects are difficult to detect. Forward and cross-decomposition models have been used to identify the causal contribution of collinear features onto multi-dimensional observations (e.g. (Naselaris et al., 2011) ). (Kriegeskorte et al., 2008) , and consists in analyzing the similarity of brain responses associated with specific categorical conditions (e.g. distinct images), by (1) fitting one-against-all classifiers on each condition and (2) testing whether these classifiers can discriminate all other conditions. The resulting confusion matrix is then analyzed in an unsupervised manner to reveal which conditions lead to similar brain activity patterns. B2B differs from RSA in that (1) it uses regressions instead of classifications, and can thus generalize to new items and new contexts and

(2) it is fully supervised.

Finally, CCA has been used in neuroimaging for a variety of purposes such as denoising and subject alignment (Hotelling, 1936; de Cheveigne et al., 2019) . While CCA relates to B2B, these two methods diverge in several ways. First, CCA and B2B have different objectives: CCA aims to find the potentially numerous and poorly interpretable components where X and Y are maximally correlated, whereas B2B aims to recover the causal factors from X to Y. Second, B2B is not symmetric between X and Y : it aims to identify specific causal features by first optimizing over the decoders G and then over H. By contrast, CCA is symmetric between X and Y , and aims to find G and H such that they project X and Y on maximally correlated dimensions. Third, CCA is based an eigen decomposition of XH and Y G -the corresponding canonical components are thus mixing the X features in way that limit interpretability and potentially dilute the impact of each feature onto multiple components.

In contrast B2B assesses each feature X j on a single Y component specifically selected to maximize signal-to-noise ratio of that feature j. Fourth, and unlike B2B, CCA does not separately optimize two distinct regularization parameters for G and H. Finally, CCA does not use different data splits to estimate G and H. Together, these differences may explain why B2B reliably outperform CCA on estimating causal influences (Figs. 2 and 5 ).

5. CONCLUSION:

In this work, we proposed Back-to-Back (B2B) regression, a linear method to measure the causal influence of a potential set of variables generating multidimensional observations. B2B performs two successive multidimensional regressions: one from the output domain, and another one from the input domain. We provided a theoretical guarantee about the consistency of B2B, and compared it to several baselines in controlled synthetic experiments. We also applied B2B to a recent brain imaging dataset, analyzing the timing of brain responses and their connection to word features. We obtained results consistent with prior work in neuroscience literature, confirming the efficacy of B2B for real data analysis. Proof. Let Ĝ and Ĥ be the solutions of the first and second regressions of B2B.

A APPENDIX:

Since Ĝ is the least square estimator of

X from Y Ĝ = arg min G E[ Y G -X 2 ]

Replacing Y by its model definition Y = (XE + N )F , we have

Ĝ = arg min G E[ X -(XE + N )F G 2 ] = arg min G E[ X -XEF G + N F G 2 ]

Since N is centered and independent of X, we have Let us now prove that EF Ĝ = F Ĝ.

Let F † be the pseudo inverse of F , and Z = F † EF Ĝ, we have F Z = F F † EF Ĝ

Since F is full rank on Img(E), we have F F † E = E, and F Z = EF Ĝ

As E is a binary diagonal matrix, it is an orthogonal projection and therefore a contraction, thus

N EF Ĝ 2 ≤ N F Ĝ 2 and X -XEF Z 2 + N F Z 2 = X -XEF Ĝ 2 + N EF Ĝ 2 ≤ X -XEF Ĝ 2 + N F Ĝ 2

But since Ĝ = arg min G X -XEF G 2 + N F G 2 , we also have To test whether each of the methods robustly scales to an increasingly large number of potential causes X, we enhanced the four ad-hoc features (word length, word frequency, word function, dummy variable) with another ten features. These additional features corresponds to the first dimensions of word embedding as provided by Spacy (Honnibal & Montani, 2017) . The results shown in Figure 7 , show that the feature importance of ad-hoc features as derived by B2B remain unchanged and are actually improved.

X -XEF Ĝ 2 + N F Ĝ 2 ≤ X -XEF Z 2 + N F Z

:



Figures and Tables:

Figure 3: Ninety subjects each read approximately 2,700 words while their brain activity was recorded with MEG.Top. Average brain response to words (word onset at t=0 ms), as viewed from above the head (red= higher gradient of magnetic flux). Bottom. Each line represents a magnetometer, color-coded by its spatial position. Posterior responses, typical of primary visual cortex activity, peak around 100 ms after word onset and are followed by an anterior propagation of activity typical of semantic processing in the associative cortices.

Figure 4: Multiple models (color-coded) are compared on their ability to reliably predict single-trial MEG signals evoked by words. Left. Average improvement of correlation coefficient ∆R for each of the four features (rows). Error bars indicate standard error of the mean (SEM) across subjects. Right. Average ∆R across time for each subject (dots). Top horizontal lines indicate when B2B significantly outperforms other methods (red) and vice versa.

-Y Ĝ 2 ] = arg min H E[ XH -(XE + N )F Ĝ 2 ] = arg min H E[ X(H -EF Ĝ) 2 ] + E[ N F Ĝ 2 ] = arg min H E[ X(H -EF Ĝ) 2 ]a positive quantity which reaches a minimum (zero) for Ĥ = EF Ĝ (15)

Figure7: Comparison of ∆R when the models are tested on four variables (top) and when the models are tested on an these four variables as well as another 10 word-embedding features (bottom). These results illustrate that, unlike Regularized CCA, B2B remains robust even when the number of tested factors increases.

Peer review (iclr2020):
The authors introduce a method (B2B) for disentangling effects of correlated predictors in the context of high dimensional outcomes. Their model assumes the outcomes are constructed by a linear transformation of a set of true causes plus measurement noise. Specifically, they fit the model Y=(XE+N)F, where X are the predictors, E is a binary matrix indicating the true causes, N is a noise term, and F is a mixing term. They provide a closed form solution the model fit based on a pair of l2-regularized regressions. They simulate from the given model and provide comparisons against least-squares regression, canonical correlation analysis, and partial least squares. They also apply the method to brain imaging data containing 102 individuals reading 120 sentences plus scrambled sentences, with the goal of inferring which features of the words have an effect on imaging results. This paper appears to be technically sound, but it should be rejected based on 1) the relatively limited applicability of the model and 2) a lack of thorough experimentation indicating that this is an appropriate method under more general circumstances. It is odd that the model assumes the outcomes are measured without error. Instead, it is assumed that the causes are measured with error, and mixed via F. A more appropriate model may be: Y=(XE+N)F+M, where an additional noise term M allows for Y to be measured imprecisely. Viewed in this light, the model starts to look a lot like canonical correlation analysis. Consider a model Y=ZF+M, X=ZG+N, if dim(Z) = dim(X) and G is invertible, this can be re-written as Y=(X inv(G)-N)A+M, and we arrive at a similar model with specific restrictions about the structure of inv(G) (E will in general not be invertible, so they are not the same). It is particularly odd, then, that the authors do no comparisons against any regularized form of CCA, which would seem to be the most natural method to use in this circumstance. Moreover, in the simulations where they show B2B outperforms CCA, they use 1000 training samples with 10-100 possible causes. In their experiments on real data, where it seems the CCA results are much closer to the results that B2B gives, they have 2700 samples and 4 possible causes. This seems to imply the real data case might be better conditioned than the simulated case, so that regularization would have less of an impact. In conclusion, the authors present a sound method for disentangling correlated possible causes when the outcome is high-dimensional. However, the authors do not provide enough evidence that this method is generally useful and better than established methods to merit acceptance to ICLR. A comparison to regularized CCA, application to more datasets and simulations under violations of the model would greatly improve the paper. I also have two minor points. 1) the word “causal” can mean many things, and here it refers specifically to disentangling correlated predictors, rather than confounding in observations or direction of effect. It would improve the paper to add some discussion of this point. 2) The comparisons in the experiments are done between E estimated from B2B and E=sum_j H_j^2 for other methods that do not directly estimate E. However, a more natural comparison might be against EF as this also includes estimates of the strength of influence of each observation, which is implicit in the sum above.