Measuring causal influence with back-to-back regression: the linear case
Authors: Jean-Remi King
Francois Charton
Maxime Oquab
David Lopez-Paz
Venue:  (2020)

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.


Under review as a conference paper at ICLR 2020
MEASURING CAUSAL INFLUENCE WITH BACK- TO-K-TO-BAC REGRESSION: THE LINEAR CASE
Anonymous authors Paper under double-blind review
ABSTRACT
Identifying causes from observations be particularly challenging when poten- tial factors are difficult to manipulate individually and ii) observations are complex and multi-dimensional. To address this issue, we introduce "Back-to-Back" regres- sion (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 iquares regression and cross- deooppicooooiiiionn 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 respec- tively, more reliably detected with B2B than with other standard techniques.
INTRODUCTION
Natural sciences are tasked to find, from a set of hypothetical factors, the minimal subset that suffices 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 word presented on a subject's 's retina modulates brain activity can be surprisingly difficult. Indeed, the frequency of words in natural language covaries with other factors such 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 set of words that control for all of these factors simul taneousl, it is thus common to use forward modeling, i.c. 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 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 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, 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
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 X and YG are maximially correlated or maximally covarying respectively. Because CCA and PLS are based on 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 al. 1995).

---

Under review as a conference paper at ICLR 2020
Observations
Factors
Cause selection
Noise
Cause-effect mapping
Y E Baxd,
X Raxd
E Dd.xd
N Ruxa,
F P..xX
X: regression from to
X P(X) ~ P(N)
É = =X xx Y2 Y
2) regression from X X
Figure 1: Back-to-back regression identifies the subset factors E X that influence some observations Y by 1) regressing from Y X to obtain X, and 2) returning the diagonal of the regression coefficients from to X.
Here, we introduce the "back-to-bace 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 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 temporal-delays and/or nonlinear interactions in systems where the directionality of causation (from to Y or vice versa) is unknown (e.g. (Peters
al. 2017
Granger 1969 Janzing al. 2013 Schölkopf al. 2016).
we show with synthetic data that it outperforms state-of-the-an art forward. backward and cross cross-dcoomposition techniques in identifying causal influence (Section B.Ip. 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
After detailing B2B method and proving convergence (Section
BACK-TO-BACK REGRESSION
We consider the measurement of multivariate signal Y E Raxd, generated from set of putative causes X E Raxds, via some unknown linear apparatus F E Rd, Not all the variables in X exert causal influence on Y. By considering square binary diagonal matrix of causal influences E E Xds, denote by X the causal factors of Y. In summary, the problem can be formalized as:
y: = (z, F
(1)
where i is given sample, and n; is sample-specific noise drawen from centered distribution. While the triplet of variables and N are independent, we allow each of them to have any form of covariance. In practice, we observe samples (X, Y) from the model. This problem space, along with the sizes all variables involved, is illustrated Figure T. Given the model Equation (T).
the goal Back-to-Back Regression (B2B) is 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 to X, and construct the predictions X Y&. This backward regression recovers the correlations between Y and each factor of X. Second, estimate the linear regression
2

---

Under review as a conference paper at ICLR 2020
coefficients H from X to X. The diagonal of the regression coefficients H. denoted by E diag H), the desired estimate of the causal influence matrix E, as detailed the A.1
Rifkin & Lippert 2007 B2B has closed form If using 12-reuularized least-squares (Hoer solution:
1959
G (Y + -y = + xxx
(2)
(3)
where Ax and Ay 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 result in overfitting, spurious correlations the data absorbed by the first regression will be leveraged by the second one. To avoid this issue, split sample (X, Y) into two splits (X1, Y1) and (X2, Y2). Then, the first regression is performed using (X Yi), and the second regression is performed using (X2, Y2)- 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 valsdvaliootion scheme as detailed in (Rifkin Lippert 2007). as follows:
= EEx - diag(Ex / diag(/ ExG)
(element-wise division)
(4)
where Ex is the X kernel matrix and where G is computed with an eigen decomposition of X:
Σx QVQ
(5)
G = Q(V + --)
where Q. and À are the eigen vectors, eigen values and regularization, respectively.
We summarize the B2B procedure in Algorithm guarantee on the correctness of B2B.
[ The rest of this section provides theoretical
Algorithm 1: Back-to-back regression.
Input: input data X [_n output data e gaxd, number of repetitions m N. Output: estimate of causal influences Ddd xd,
c + 0;
1
for i .......2 do (X, Y) + ShufflcRows((X, Y)):
3
(X,, Y1), (X2, + SpliRRowinHHal((X
5
LinearRegression (Y1, x
G = (YTY, + -yy
A LigearRegression Y2G):
Ā = (X2 X2 +Ax) -x, Y2G
E t E + diag( ff;;
7
8 end
E + B/m:
W LinearRegression(XX
return É, W
2.2 THEORETICAL GUARANTEES
Theorem (B2B consistency general case). Consider the B2B model from Equation Y = (XE +
N)F, N centred and full rank noise. Let Img(Mr refers to the image of the matrix M. If F and
X are full-rank on the Img(E), then, the solution of B2B, H, will minimize ming ||X XH1
||N H1 satisfy E.f. H
Proof. See Appendix A.I
3

---

Under review as a conference paper at ICLR 2020
Since EH- = H, we have
= arg min H
X XEHII ||NEI| (EX XE+ EN NE)EXXX
(6)
Assuming, without loss of generality, that the active features in E are the k & k [0, dg) first features, and rewriting X = (X1. X2) and N = (N1.NN) (X1 and Ni containing the ls first features). we have:
Ex,x, Ex,X2 Ex,x2
NN,N ENINs NN.
)
x
N
Xx,X2
INaNs
where AAB is the covariance of A and B, and:
i = 1=( ((Exxx + -.NN)) (haNn))''An Exxx, (Xxxxx + NN, -
(8)
diagı (fr) diag((Ex,,, + ENGN) -xxx,) diag(( Exix, ENINN, In the absence of noise, have EN,N, and diag, (ft) and diag( H diag(F) F.)
(9)
Therefore, we recover E from H.
In the presence of noise, the causal factors of E correspond to the positive elements of diag(ff), (ft). The methods recover them are presented in the Appendix A.4
EXPERIMENTS
perform two sets of experiments to evaluate B2B: one on controlled synthetic data, and a second one real, large-scale maggett--nncehaalggrabh (MEG) dataset. We use scikit-learn's 's PLS and RidgeCV (Pedregosa et al., 2011 as well as pyreca'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 to evaluate the ability of B2B terms of prediction of independent and identically distributed data, as well method 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 scalar that modulates the signal-to-noise ratio. Here, F E Xdy contains entries drawn from A A(O0, o²) where o² inversely proportional to de, X E Raxd, contains rows drawn from (0, Σx), N E [xdd contains rows drawn from (0, Ey), E E Xd= a binary diagonal matrix containing Ne ones, Σx AAT where A Rd xd contains entries drawn from (000 ΣN BB where B E Rdd xd, contains entries drawn from A (0, o²), and the factor E R...
To simulate wide range of experimental conditions, we sample 10 values in log-space for de, dy [10, 100). n. E [3, 63). h E [0.001, 10). We discard the cases where n. > dg. limit dg dy to 100 to keep the running time under 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 al. 2011 and pyreca (Bilenko & Gallant 2016):
3.2 BASELINE MODELS
Forward regression consists of an 12-regularized "ridge" regression from the putative causes X to the observations Y:
Hywd (xT X -XX-
(10)
4

---

Under review as a conference paper at ICLR 2020
1.0
Forward
RegCCA
PLS
B2B
AUC
0.8
0.6
10
10
10a
10a
10
10a
101
10a
snr
ne
d2
d
Forward
RegCCA
PLS
B2B
0.4
0.2
10-3
10
10a
101
10a
shr
dx
Figure 2: Synthetic experiments. Average AUC (top) and Feature Importance AR (bottom) when varying experimental conditions individually. Higher is better. B2B compares favorably in all cases.
Backward regression consists of an 12-regularized "ridge" regression from Y to X:
Gowd = T + xXx
(11)
CCA finds Geca E Rddadd and Heca E Rdddda s.t. X and are maximally correlated a latent Z space:
Geeas Heea arrmax YYG)
(12)
G.H
PLS finds Gpts E Rd and Hpts Rdaide s.t. X and Y are maximally covarying in a latent Z space:
Gyte: Hpts = argmax ( YG")
(13)
We employ 5-fold cross-validation to select the optimal number of components for CCA and PL.S. Regressions were /2-regularized with À regularization parameters fitted with the efficient leave-one- out procedure implemented in scikit-learn RidgeCV (Pedregosa et al. 201T.
3.3 EVALUATING CAUSAL DISCOVERY FROM MODELS' COFFFICIENTS
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 loading vector H, per feature i (or vector G for the backward model). To transform such vector into estimated causal contribution É, take the sum of square coefficients: É. E, H!
To estimate whether models accurately identify causal factors, we compute the area-under-the-curvee (AUC) across factors AUC É). The AUC allows evaluating the capacity of models at detecting the causal importance of factors when ground truth labels are available, the case in this setup.
We report AUC results Figures (top) and (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 thus not be estimated. To address this issue, 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.e. "knock-out X'). This procedure

---

Under review as a conference paper at ICLR 2020
350 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 t=0 ms), as viewed from above the head (red= higher gradient of magnetic flux). Bottom. Each line rep- resents 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 anterior propagation of activity typical of semantic processing in the associative cortices.
Nave=90 90
100
50
-100
500
600
Time from word onset (ms)
results in two correlation metrics Ryall and Rknockout- whose difference AR, = Rfull - Rknockout indicates how much each X, improves the prediction of Y. In figures, AR is the average of A R. A higher score means that for prediction, the model relies on individual features rather than combinations of features.
We show in Appendix A.33 pseudo-code to assess feature importance for our algorithm as well as baselines. For the Backward Model, feature importance cannot be assessed the X colinearity is never taken into account.
We show in Figures (bottom) and (right, in Appendix) that method outperforms baselines.
3.5 MAGNETORNCEPHALOCGROO
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 (MEG) recordings of 102 healthy native-Dutch adults who participated in reading task. Twelve subjects were excluded from the analysis because of corrupted file headers. Subjects were exposed 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 351 ms average (min: 300ms, max: 1400ms). Successive words were separated by 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 parame- ters (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. H. and 10Hz for the low and high cut-off frequencies. The raw data was then segmented 100ms before word onset and Is after word onset (t Oms corresponds to word onset). Finally, each resulting segment was baseline-corrected between 100ms and Oms, 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 E Raxdd of n Rs 2700 words by dy = 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 brain responses. We consider four distinct but colinear features. First, 'Word Length' refers to the total number of letters. Word Length expected to specifically cause a variation in the early evoked MEG responses (i.e. from 100 ms after stimulus onset) elicited by the retinotopicalll-tum 2.14).). 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
visual cortices (e.g. (Pegado al.

---

Under review as a conference paper at ICLR 2020
package (Speer et al. 2018 Word Frequency is expected to specifically cause variation in the late evoked MEG responses (i.e. from 400 ms), because 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, conjunction, a determinant, pronoun or 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 dummy feature, constructed from a noisy combination of Word Length and Word Frequency: durany = z(length) + = ('requency) + N, where normalizes features and A is random vector sampling Gaussian distribution (all terms thus have a zero-mean and a unit-variance) This procedure yields an X E Raxdz matrix of n ss 2700 words by d 4 features for each subject. Each of the columns of is normalized to have a mean and a standard deviation of and respectively.
3.5.3 MODELS AND STATISTICS
We compare B2B to four standard methods: Forward regression, Backward regression, CCA and PL.S, as implemented in scikit-learn (Pedregosa et al., 2011) and (Bilenko & Gallant 2016), and optimized with nested voss-validation over twenty 12 regularization parameters logarithmically spaced between 10-4 and 10t (for regression and CCA methods) or to 4 canonical components (for PLS).
We used the feature importance described in Algorithm 2 to assess the extent to which each feature X, 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 performed using two-sided Wilcoxon test across subjects (n=90) using the average AR R across time. Corresponding effect sizes are shown in Figure and p-values 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 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 A (i.e. the improvement of Y prediction induced by the introduction of 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 AR < 0 (all p .089).
Figure 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: AR for all models (all < 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. 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 < 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.
7

---

Under review as a conference paper at ICLR 2020
4
RELATED WORK
Forward PLS
RegCCA
B2B
Word Length
0.15
Forward and cross-decomposition models have been used to identify the causal contribution of collinear fea- tures onto multi-dimensional observa- tions (e.g. (Naselaris et al. 2011). These approaches typically lead to multiple coefficients for each features (i.e. one per dimension of or one component respectively). Further- more, these coefficients can be diffi- cult to summarize into single causal estimate. By contrast, B2B quickly (Fig. 6 leads to single unbiased scalar values E tending towards and O for causal and non-causal features respectively.
AR
AR
0.00
0.00
0
500
1000
Word Frequency
0.06
AR
AR
0.00
0.00
0
500
1000
Word Function
0.06
A variety of other statistical methods applied to neuroimaging data have been proposed to clarify what is be- ing represented in brain responses Le. what feature causes specific brain activity. One of the popular linear method is Representational Similarity Analysis (RSA) (Kriegeskorte et al., 2008), and consists in analyzing the similarity brain responses associ- ated with specific categorical condi- tions (e.g. distinct images), by (I) fitting one-against-all classifiers on each condition and (2) testing whether these classifiers can discriminate all other conditions. The resulting con- fusion matrix is then analyzed in an unsupervised manner to reveal which conditions lead to similar brain activ- ity patterns. B2B differs from RSA in that (1) it uses regressions instead of classifications, and can thus general- ize to new items and new contexts and (2) it is fully supervised.
AR
AR
0.00
500
1000
Dummy Variable
0.06
AR
0.00
0.00
0
500
1000
Model
Time (ms)
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 corre- lation coefficient AR for each of the four features (rows). Error bars indicate standard error of the mean (SEM) across subjects. Right. Average AR across time for each subject (dots). Top horizontal lines indicate when B2B significantly outperforms other methods (red) and vice versa.
Finally, CCA has been used in neuroimaging for a variety of purposes such as denoising and subject alignment (Hotelling 1936 de Cheveigne 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 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 H and YG 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 on single Y component specifically selected to maximize signal-to-noise 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. and 5)
8

---

Under review as a conference paper at ICLR 2020
5 CONCLUSION
In this work, we proposed Back-to-Back (B2B) regression, linear method to measure the causal influence of potential 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 theoretical guarantee about the consistency of B2B, and compared it to several baselines in controlled synthetic experiments. We also applied B2B to 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.

---

Under review as a conference paper at ICLR 2020
REFERENCES
Natalia Y Bilenko and Jack I. Gallant. Pyreea: regularized kernel canonical correlation analysis in python and applications to neuroimaging. Frontiers in neuroinformattis, 10:49, 2016.
Leo Breiman. Bagging predictors. Machine learning, 24(2):123-140, 1996.
Radoslaw Martin Cichy, Dimitrios Pantazis, and Aude Oliva. Resolving human object recognition in space and time. Nature neuroscience, 17(3):455, 2014.
Alain de Cheveigne, Giovanni M Di Liberto, Dorothee Arzounian, Daniel DE Wong, Jens Hjortkjaer, Søren Fuglsang, and Lucas C Parra. Multiway canonical correlation analysis of brain data. Neurolmage, 186:728-740. 2019.
Karl J Friston, Andrew P Holmes, Keith J Worsley, J-P Poline, Chris D Frith, and Richard SJ Frackowiak. Statistical parametric maps in functional imaging: general linear approach. Human brain mapping, 2(4):189-210, 1994.
Alexandre Gramfort, Martin Luessi, Eric Larson, Denis A Engemann, Daniel Strohmeier, Christian Brodbeck, Roman Goj, Mainak Jas, Teon Brooks, Lauri Parkkonen, et al. Meg and eeg data analysis with mne-python. Frontiers in neuroscience, 7:267, 2013.
Alexandre Gramfort, Martin Luessi, Erie Larson, Denis A Engemann, Daniel Strohmeier, Christian Brodbeck, Lauri Parkkonen, and Matti S Hämäläinen. Mne software for processing meg and eeg data. Neuroimage, 86:446-460, 2014.
Clive WJ Granger. Investigating causal relations by econometric models and cross-spectral methods. Econometrica: Journal of the Econometric Society, pp. 424 438, 1969.
Martin N Hebart and Chris I Baker. Deconstructing multivariate decoding for the study of brain function. Neuroimage, 180:4-88, 2018.
Arthur E Hoerl. Optimum solution of many variables equations. Chemical Engineering Progress, 55 (11):69-78, 1959.
Matthew Honnibal and Ines Montani. spacy 2: Natural language understanding with bloom embed- dings, convolutional neural networks and incremental parsing. To appear, 2017.
H. Hotelling. Relations between two sets of variables. Biometrika, 129):129-149, 1936.
Alexander Huth, Wendy A de Heer, Thomas L Griffiths, Frédéric E Theunissen, and Jack L Gallant. Natural speech reveals the semantic maps that tile human cerebral cortex. Nature, 532(7600):453, 2016.
Dominik Janzing, David Balduzzi, Moritz Grosse Grosse-Wentruy Wentrup, Bernhard Schölkopf, et al. Quantifying causal influences. The Annals of Statisties, 41(5):2324-2358, 2013.
G. V. Kass. Significance testing in automatic interaction detection (a.i.d.). Journal of the Royal Statistical Society. Series C (Applied Statistics), 1482):178-88, 1975.
Jean-Rémi King, Laura Gwilliams, Chris Holdgraf, Jona Sassenhagen, Alexandre Barachant, Denis Engemann, Eric Larson, and Alexandre Gramfort. Encoding and decoding neuronal dynamics: Methodological framework to uncover the algorithms of cognition. 2018.
Nikolaus Kriegeskorte, Marieke Mur, and Peter A Bandettini. Representational similarity analysis- connecting the branches of systems neuroscience. Frontiers in systems neuroscience, 2:4, 2008.
Marta Kutas and Kara D Federmeier. Thirty years and counting: finding meaning in the n400 component of the event-related brain potential (erp). Annual review of psychology, 62:621-647, 2011.
Ludovic Lebart, Alain Morineau, and Marie Piron. Statistique exploratoire multidimennsommelle. volume 3. Dunod Paris, 1995.
10

---

Under review as a conference paper at ICLR 2020
J. A. Morgan and J. N. Sonquist. Problems in the analysis of survey data: and proposal. J. Amer. Statist. Ass., (58):415 1963.
Thomas Naselaris, Kendrick N Kay, Shinji Nishimoto, and Jack L Gallant. Encoding and decoding in fmri. Neuroimage, 56(2):400410, 2011.
Kenneth Norman, Sean M Polyn, Greg Detre, and James V Haxby. Beyond mind-reading: voxti pattern analysis of fmri data. Trends in cognitive sciences, 10(9))424 430, 2006.
Fabian Pedregosa, Gačl Varoquaux. Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathicu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al. Scikit-learn Machine learning in python. Journal of machine learning research, 12(0ct):282--2830, 2011.
Felipe Pegado, Enio Comerlato, Fabricio Ventura, Antoinette Jobert, Kimihiro Nakamura, Marco Buiatti, Paulo Ventura, Ghislaine Dehaene-Lambertz, Régine Kolinsky, José Morais, et al. Timing the impact of literacy on visual processing. Proceedings of the National Academy of Sciences, 111 (49): 2014.
Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of causal inference: foundations and learning algorithms. MIT press, 2017.
Ryan M Rifkin and Ross A Lippert. Notes on regularized least squares. Technical report, MIT, 2007.
Jan-Mathijs Schoffelen, Robert Oostenveld, Nietzsche HL Lam, Julia Uddén, Annika Hultén, and Peter Hagoort. A 204-subject multimodal neuroimaging dataset to study language processing. Sciietific data, 671):17, 2019.
Bernhard Schölkopf, David W Hogg, Dun Wang, Daniel Foreman-Mackey. Dominik Janzing, Carl- Johann Simon-Gabriel, and Jonas Peters. Modeling confounding by half-sibling regression. Proceedings of the National Academy of Sciences, 113(27):7311-7398 2016.
Robyn Speer, Joshua Chin, Andrew Lin, Sara Jewett, and Lance Nathan. Luminnsoinsighrwooddnng v2.2, October 2018. URL httlthttss://cco wtt.//o...ogg/../o/ 5281/ 2tt////eeooo.1111 14 1443582
Walter JB Van Heuven, Pawel Mandera, Emmanuel Keuleers, and Mare Brysbuert. Subtlex-ut -kklee--uk: A new and improved word frequeney database for british english. The Quarterly Journal of Experimental Psychology, 67(6): 76-1190, 2014.
Sebastian Weichwald, Timm Meyer, Ozan Özdenizci, Bernhard Schölkopf, Tonio Ball, and Moritz Grosse Grosse-Wentrup. Causal interpretation rules for encoding and decoding models in neuroimaging. Neurolmage, 110:48-50 2015.

---

Under review as a conference paper at ICLR 2020
A
APPENDIX
A. PROOF OF CONSISTENCY THEOREM
Proof of the theorem in 2.2
Theorem 2 (B2B consistency general case). Consider the B2B model from equation
+
with N centred and full rank noise.
If F and X are full-rank on Img(E), then, the solution of B2B, H minimizes
min - - XH N
and satisfies
En = A
Proof. Let G and Á be the solutions of the first and second regressions B2B.
Since G is the least square estimator X from Y
-= arg min B[| x G
Replacing Y by its model definition Y (XE + N)F, we have
= arg min E|| (XE - FNG arg min B|J XEFG+
Since N is centered and independent of X, we have
G = arg min
| XEFGI| |NFG]]
(14)
G
Samely, for H, have
A arg min E[]| ------YYII Y
arg min B| (XE N)F
arg min B X(H - -
arg min B X(H 1F
#
a positive quantity which reaches a minimum (zero) for
=
(15)
Let us now prove that EFG = FG.
Let Ft be the pseudo inverse of F, and Z = FFEFG, we have FZ = FFFEFG
Since F' is full rank on Img(E), we have FF = E, and FZ = EFG
As E is binary diagonal matrix, is an orthogonal projection and therefore contraction, thus
|NEFG < ||NF
and ||X
||X - -IX
But since & arg ming
|| -EFGGI |INFG| also have
|x x (
< X - XEFZ
12

---

Under review as a conference paper at ICLR 2020
Summarizing the above,
x - T-XFFo| + | XEF + & ||X I
X xx-x
NFG
|| - INEFGI
NNG |NEFG
N being full rank, this yields = Replacing into (14) and setting H = F.G,, we have
G arg min X KEFGI'+INNGII
arg min
X XEFGII+IINERR
A = arg min |X
x |
Finally, EA = EEF = EFC- = H, since E, a binary diagonal matrix, involutive. This completes the proof.
13

---

Under review as a conference paper at ICLR 2020
A.2 MODELING MEASUREMENT NOISE
Equation does not explicitly contain a measurement noise term. Yet, in moast experimental cases, the problem is best described as:
Y= (XE+ NNF++ M
(16)
with M E Rexd,
This equation actually equivalent to Equation given our hypotheses. Indeed, we can rewrite
M = MF-IF over Img(F) F), which leads to:
Y (XE+ M NXE+ FF -
Consequently, assuming that F full rank on Img(X (X F), B2B yields the same solutions to equations
andd
A.3 FEATURE IMPORTANCE
For B2B, feature importance is assessed as follows:
Algorithm 2: B2B feature importance.
Input: X traann E Raxd, Xtest R=''dd, Ytrain E Raxdy Yeest Rr'xd, Output: estimate of prediction improvement AR E DMd.
1 H, G B2B (troin B2B(Xtroin: Yirain);
2 Rfull = corr( Xtest H, G);
3 for 1,..., dr do
K Id;
5
[[i] f- 0:
R2 corr(Xieed
Xtest K KH, Ytest
test
7
AR, = Riull Rai
8 end
return AR
For the Forward Model, the feature importance is assessed as follows:
Algorithm 3: Forward feature importance.
Input: Xtrain E Raxdd, Xtest E TR' xdz Ytrain E RRaxd, Ytest E RR xdy Output: estimate of prediction improvement AR Dydzidy
1 H = AenearrRgression(X..... (Xtrinn Ytrain) Ryall corr(X Yrest);
2 for 1,..., de do
K = Id:
K[i] + 0;
Rg corT( X.r(X KH, Ytest);
AR, R&guI Re:
7 end
8 return AR
14

---

Under review as a conference paper at ICLR 2020
For the CCA and PLS models, the feature importance is assessed follows:
Algorithm 4: CCA and PLS feature importance.
Input: Xtrain E Raxd, Xtest R"''dd, xd,, Ytrain E Raxd, Yest E Rr'xdd, Output: estimate of prediction improvement AR E Ddzids.
1 H, G CCACCCAAAAA traint Virain);
2 Rfull corr( (ten(teetlll G);
3 for 1....,dy do
K Id:
5
K[i] f- 0:
Rg corT( X test K KH, G);
7
AR; Ri Re:
8 end
return AR
For the Backward Model, feature importance cannot be assessed because there no prediction.
A.4 RECOVERING
In case noise, B2B yields non binary E. Three thresholding rules can be used binarize its values thus explicitly recover causal features.
First, given known signal-to-nn -ggnae-t-nooie ratio,the threshold above which a feature should considered be causal can be derived analytically. Indeed, Equation 9 implies that the & first diagonal elements of H are bounded:
OX S diag, (f) Ox, + ON,
OX
S
ax, + +ON
where x,, xx., ON, and d,, denote the largest and smallest eigenvalues of X, and ΣN,
The average value μ of non-zero coefficients diag(f) is the trace of Á divided by k and can be computed as
Var(X) Var(X) Var(N)
(17)
The decision threshold between causal and non-causal elements is thus a fraction Jr, whose proportion arbitrarily depends on the necessity to favor type and type errors. In practice, we cannot use this procedure for our MEG study, because signal-to-noise ratio is unknown.
Second, diag( ( can be binarized with the Nonquist-Morgant criterion (Morgan & Sonquist 1963). ) non-aarametric clustering procedure separating small and large values in given set. This procedure maximizes the ratio of inter-group variance while minimizing the intra-group variance, over all possible splits of the diagonal into p largest values and d, p smallest values. Let ma and m, be the average values of the two clusters, and dz p their size, and the total variance of the sample, Sonquist-Morgan criterion maximizes Kass 1975p:
p(d, p) (m1 m2)
(18)
dg
This procedure assumes that there exists at least one causal and at least one non-causal feature. Third, second-order statistics across multiple datasets can be used to identify the elements of diag (it that sgnificantly different from O. This procedure is detailed in the method section of our MEG experiment.
Overall, these three procedure thus varies in their additional assumptions: i.c. (1) known signal-to- noise ratio, (2) the existence of both causal and non-causal factors or (3) independent repetitions of the experiment.
15

---

Under review as a conference paper at ICLR 2020
AUC
AUC
AUC
1.0
1.0
1.0
0.8
0.8
0.8
B2B
B2B
B2B
0.6
0.6
0.6
0.4
0.4
0.4
0.5
1.0
0.5
1.0
0.5
1.0
Forward
PLS
RegCCA
AR
AR
AR
0.04
0.04
0.04
B2B
0.02
B2B
0.02
B2B
0.02
0.00
0.00
0.00
0.000
0.025
0.000
0.025 PLS
0.000
0.025 RegCCA
Forward
Figure Synthetic experiments. Distribution (over conditions) AUC (top) and Feature Importance AR (bottom) metrics between method (y-axis) and the baselines (x-axis). Each dot is distinct synthetic experiment. Dots below the diagonal indicates that B2B outperform the tested model.
wall-clock time (seconds)
10
100
10-1
10-2
Forward PLS RegCCA B2B
Figure 6: Wall-clock run-time for our method B2B and for the baselines. Each dot is a distinct synthetic experiment. B2B runs much faster than cross- decompsccoooiiitioon baselines.
16

---

Under review as a conference paper at ICLR 2020
Word Length
Word Frequency
model
Word Function
PLS
RegECA
Dummy Variable
828
0.01
0.04
Length
Word Frequency
Word Function
Dummy
embedding
embedding
effect
embedding
embedding
embedding
embedding
embedding
embedding
model
Forward
embedding
4,
RegCCA
embedding
828
0.00
0.01
0.02
0.03
0.04
Figure 7: Comparison AR when the models are tested on four variables (top) and when the models tested an these four variables 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.
B ADDITIONAL FIGURES
B. ROBUSTNESS TO INCREASING NUMBER OF FACTORS
To test whether each of the methods robustly scales increasingly large number of potential causes X, we enhanced the four ad-hoc features (word length, word frequency, word function, dummy variable) with another features. These additional features corresponds 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 actually improved.
17


Reviews:
[{"review_id": "rJl4BEwAYH", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "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.\n\nThis 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.\n\nIn 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.\n", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}, {"review_id": "HyxooSwIYB", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "\nThe paper provides an iterative linear method to identify causal influences of putative cause matrix to signal matrix. The idea is a natural extension of previous forward and backward such as CCA and PLS. The paper has provided consistency guarantee and several synthetic and real data experiments as support. \n\nTechnical questions:\n\n1) The estimation of binary causal influence matrix E is set as \\hat{E} = Diag(\\hat{H}). Why is Diag(\\hat{H}) guaranteed to have binary elements?\n\n2) The Theorem 1 needs more explanation about why it proves consistency, which is currently isolated from other parts of the paper. Why E\\hat{H} = \\hat{H} guarantees the consistency of \\hat{E} = Diag(\\hat{H})?  For example, \\hat{H} can have more all-zero rows than E which still satisfies E\\hat{H} = \\hat{H}  but  \\hat{E}  is not equals to E. In an extreme case,  \\hat{H} = 0 will have  E\\hat{H} = \\hat{H} but is clearly not consistent. \n\n3) How does the Eq (4),(5) give an estimation \\hat{E}?\n\n4) By Eq. (10,11), H and G seem to be determined give X,Y. Then what are the maximization of Eq. (12,13) over? \n\n\nGeneral comments:\n\n1) How does the problem in Eq.(1) differ from variable selection in linear regression where a plenty algorithms exist such as LASSO, spike-slab prior, SCAD, etc. ?\n\n2) The experiments are a bit weak with a simple synthetic experiment and a real dataset with just four features. Can the experiment directly demonstrate the correctness of the theorem?\n\nTypo:\n\nPage 1, “are are based on”\n\nIn general, I think B2B is an algorithm that has improvement over CCA and PLS. I am looking forward to the author response to address my above concerns.\n\n##############\nI have read the author's feedback which addresses some of the confusion parts in the paper. I maintain the current rating mainly because of the experimental strength. ", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_checking_correctness_of_experiments": "I carefully checked the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}, {"review_id": "Bkgo2vCzcr", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper proposes \"Back-to-Back\" regression for estimating the causal influence between X and Y in the linear model Y=(XE+N)F, where the E denotes a diagonal matrix of causal influences. Furthermore, this work theoretically shows the consistency of B2B and the experiments also verify the effectiveness of this approach.\n\nThe writing is well and clear and there are some minors issues:\n- Further analysis and explanation for using \\hat{E}=Diag(H)to estimate the causal influence might be needed.\n- The model defined in Fig. 1 seems the influence E should have a sparse diagonal vector. It is possible to introduce an L1 regulation in E?\n\n##############\nAfter reading the author's feedback and the comments from other reviewers, I keep the current rating but tend to a borderline score and it is ok if it must be rejected because of the concerns of limited applicability and the experimental.\n", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The authors introduce a method for disentangling effects of correlated predictors in the context of high dimensional outcomes. While the paper contains interesting ideas and has been substantially improved from its original form, the paper still does not meet the quality bar of ICLR due to its limitations in terms of limited applicability and experiments. The paper will benefit from a revision and resubmission to another venue.
Decision: Reject