Interpretability of Multivariate Brain Maps in Linear Brain Decoding: Definition, and Heuristic Quantification in Multivariate Analysis of MEG Time-Locked Effects

Brain decoding is a popular multivariate approach for hypothesis testing in neuroimaging. Linear classifiers are widely employed in the brain decoding paradigm to discriminate among experimental conditions. Then, the derived linear weights are visualized in the form of multivariate brain maps to further study spatio-temporal patterns of underlying neural activities. It is well known that the brain maps derived from weights of linear classifiers are hard to interpret because of high correlations between predictors, low signal to noise ratios, and the high dimensionality of neuroimaging data. Therefore, improving the interpretability of brain decoding approaches is of primary interest in many neuroimaging studies. Despite extensive studies of this type, at present, there is no formal definition for interpretability of multivariate brain maps. As a consequence, there is no quantitative measure for evaluating the interpretability of different brain decoding methods. In this paper, first, we present a theoretical definition of interpretability in brain decoding; we show that the interpretability of multivariate brain maps can be decomposed into their reproducibility and representativeness. Second, as an application of the proposed definition, we exemplify a heuristic for approximating the interpretability in multivariate analysis of evoked magnetoencephalography (MEG) responses. Third, we propose to combine the approximated interpretability and the generalization performance of the brain decoding into a new multi-objective criterion for model selection. Our results, for the simulated and real MEG data, show that optimizing the hyper-parameters of the regularized linear classifier based on the proposed criterion results in more informative multivariate brain maps. More importantly, the presented definition provides the theoretical background for quantitative evaluation of interpretability, and hence, facilitates the development of more effective brain decoding algorithms in the future.


INTRODUCTION
Understanding the mechanisms of the brain has been a crucial topic throughout the history of science. Ancient Greek philosophers envisaged different functionalities for the brain ranging from cooling the body to acting as the seat of the rational soul and the center of sensation (Crivellato and Ribatti, 2007). Modern cognitive science, emerging in the twentieth century, provides better insight into the brain's functionality. In cognitive science, researchers usually analyze recorded brain activity and behavioral parameters to discover the answers of where, when, and how a brain region participates in a particular cognitive process.
To answer the key questions in cognitive science, scientists often employ mass-univariate hypothesis testing methods to test scientific hypotheses on a large set of independent variables (Groppe et al., 2011a;Maris, 2012). Mass-univariate hypothesis testing is based on performing multiple tests, e.g., t-tests, one for each unit of the neuroimaging data, i.e., independent variables. The high spatial and temporal granularity of the univariate tests provides fair level of interpretability. On the down side, the high dimensionality of neuroimaging data requires a large number of tests that reduces the sensitivity of these methods after multiple comparison correction (Bzdok et al., 2016). Although techniques such as the non-parametric cluster-based permutation test (Bullmore et al., 1996;Maris and Oostenveld, 2007), by weak rather strong control of familywise error rate, offer more sensitivity, they still experience low sensitivity to brain activities that are narrowly distributed in time and space (Groppe et al., 2011a,b). The multivariate counterpart of mass-univariate analysis, known generally as multivariate pattern analysis, have the potential to overcome these deficits. Multivariate approaches are capable of identifying complex spatio-temporal interactions between different brain areas with higher sensitivity and specificity than univariate analysis (van Gerven et al., 2009), especially in group analysis of neuroimaging data (Davis et al., 2014).
Brain decoding (Haynes and Rees, 2006) is a multivariate technique that delivers a model to predict the mental state of a human subject based on the recorded brain signal. There are two potential applications for brain decoding: (1) brain-computer interfaces (BCIs) (Wolpaw et al., 2002), and (2) multivariate hypothesis testing (Bzdok et al., 2016). In the first case, a brain decoder with maximum prediction power is desired. In the second case, in addition to the prediction power, extra information on the spatio-temporal nature of a cognitive process is desired. In this study, we are interested in the second application of brain decoding that can be considered a multivariate alternative for mass-univariate hypothesis testing. Further, we mainly focus on the linear brain decoding because of its wider usage in analyzing inherently small sample size and high dimensional neuroimaging data, compared to the complex (Cox and Savoy, 2003;LaConte et al., 2005) and non-transparent (Lipton et al., 2016) non-linear models.
In Brain mapping (Kriegeskorte et al., 2006) the pre-computed quantities, e.g., univariate statistics or weights of a linear classifier, are assigned to the spatio-temporal representation of neuroimaging data in order to reveal functionally specialized brain regions which are activated by a certain cognitive task. In its multivariate form, brain mapping uses the learned parameters from brain decoding to produce brain maps, in which the engagement of different brain areas in a cognitive task is visualized. Intuitively, the interpretability of a brain decoder refers to the level of information that can be reliably derived by an expert from the resulting maps. From the cognitive neuroscience perspective, a brain map is considered interpretable if it enables a scientist to find out the answers of three key questions: "where, when, and how does a brain region contribute to a cognitive function?" In fact, a classifier only answers the question of what is the most likely label of a given unseen sample (Baehrens et al., 2010). This fact is generally known as knowledge extraction gap (Vellido et al., 2012) in the machine learning context. Thus far, many efforts have been devoted to filling the knowledge extraction gap of linear and non-linear data modeling methods in different areas such as computer vision , signal processing (Montavon et al., 2013), chemometrics (Yu et al., 2015), bioinformatics (Hansen et al., 2011), and neuroinformatics (Haufe et al., 2013). In the context of neuroimaging, the knowledge extraction gap in classification is generally known as the interpretation problem (Sabuncu, 2014;Haynes, 2015;Naselaris and Kay, 2015). Therefore, improving the interpretability of linear brain decoding and associated brain maps is a primary goal in the brain imaging literature (Strother et al., 2014). The lack of interpretability of multivariate brain maps is a direct consequence of low signal-to-noise ratios (SNRs), high dimensionality of whole-scalp recordings, high correlations among different dimensions of data, and cross-subject variability (Besserve et al., 2007;Anderson et al., 2011;Blankertz et al., 2011;Brodersen et al., 2011;Langs et al., 2011;Lemm et al., 2011;Varoquaux et al., 2012;Kauppi et al., 2013;Haufe et al., 2014a;Olivetti et al., 2014;Taulu et al., 2014;Varoquaux and Thirion, 2014;Haynes, 2015;Wang et al., 2015). At present, two main approaches are proposed to enhance the interpretability of multivariate brain maps: (1) introducing new metrics into the model selection procedure, and (2) introducing new hybrid penalty terms for regularization.
The first approach to improving the interpretability of brain decoding concentrates on the model selection procedure. Model selection is a procedure in which the best values for the hyperparameters of a model are determined . The selection process is generally performed by considering the generalization performance, i.e., the accuracy, of a model as the decisive criterion. Rasmussen et al. (2012) showed that there is a trade-off between the spatial reproducibility and the prediction accuracy of a classifier; therefore, the reliability of maps cannot be assessed merely by focusing on their prediction accuracy. To utilize this finding, they incorporated the spatial reproducibility of brain maps in the model selection procedure. An analogous approach, using a different definition of spatial reproducibility, is proposed by Conroy et al. (2013). Beside spatial reproducibility, the stability of the classifiers (Bousquet and Elisseeff, 2002) is another criterion that is used in combination with generalization performance to enhance the interpretability. For example Yu (2013) and Lim and Yu (2016) showed that incorporating the stability of models into cross-validation improves the interpretability of the estimated parameters (by linear models).
The second approach to improving the interpretability of brain decoding focuses on the underlying mechanism of regularization. The main idea behind this approach is twofold: 1) customizing the regularization terms to address the illposed nature of brain decoding problems (where the number of samples is much less than the number of features; Mørch et al., 1997;Varoquaux and Thirion, 2014), and (2) combining the structural and functional prior knowledge with the decoding process so as to enhance the neurophysiological plausibility of the models. Group Lasso (Yuan and Lin, 2006) and totalvariation penalty (Tibshirani et al., 2005) are two effective methods using this technique (Rish et al., 2014;Xing et al., 2014). Sparse penalized discriminant analysis (Grosenick et al., 2008), group-wise regularization (van Gerven et al., 2009), smoothedsparse logistic regression (de Brecht and Yamagishi, 2012), totalvariation ℓ 1 penalization (Michel et al., 2011;Gramfort et al., 2013), the graph-constrained elastic-net (Grosenick et al., 2009(Grosenick et al., , 2013, and social-sparsity  are examples of brain decoding methods in which regularization techniques are employed to improve the interpretability of linear brain decoding models.
Recently, taking a new approach to the problem, Haufe et al. questioned the interpretability of weights of linear classifiers because of the contribution of noise in the decoding process (Bießmann et al., 2012;Haufe et al., 2013Haufe et al., , 2014b. To address this problem, they proposed a procedure to convert the linear brain decoding models into their equivalent generative models. Their experiments on the simulated and fMRI/EEG data illustrate that, whereas the direct interpretation of classifier weights may cause severe misunderstanding regarding the actual underlying effect, their proposed transformation effectively provides interpretable maps. Despite the theoretical soundness, the intricate challenge of estimating the empirical covariance matrix of the small sample size neuroimaging data  limits the practical application of this method. In spite of the aforementioned efforts to improve the interpretability of brain decoding, there is still no formal definition for the interpretability of brain decoding in the literature. Therefore, the interpretability of different brain decoding methods are evaluated either qualitatively or indirectly (i.e., by means of an intermediate property). In qualitative evaluation, to show the superiority of one decoding method over the other (or a univariate map), the corresponding brain maps are compared visually in terms of smoothness, sparseness, and coherency using already known facts (see for example, Varoquaux et al., 2012). In the second approach, important factors in interpretability such as spatiotemporal reproducibility are evaluated to indirectly assess the interpretability of results (see for example, Langs et al., 2011;Rasmussen et al., 2012;Conroy et al., 2013;Kia et al., 2016). Despite partial effectiveness, there is no general consensus regarding the quantification of these intermediate criteria. For example, in the case of spatial reproducibility, different methods such as correlation (Rasmussen et al., 2012;Kia et al., 2016), dice score (Langs et al., 2011), or parameter variability (Conroy et al., 2013;Haufe et al., 2013) are used for quantifying the stability of brain maps, each of which considers different aspects of local or global reproducibility.
With the aim of filling this gap, our contribution is three-fold: (1) Assuming that the true solution of brain decoding is available, we present a theoretical definition of the interpretability. The presented definition is simply based on cosine proximity in the parameter space. Furthermore, we show that the interpretability can be decomposed into the reproducibility and the representativeness of brain maps. (2) As a proof of the concept, we exemplify a practical heuristic based on eventrelated fields for quantifying the interpretability of brain maps in time-locked analysis of MEG data. (3) Finally, we propose the combination of the interpretability and the performance of the brain decoding as a new Pareto optimal multi-objective criterion for model selection. We experimentally, on both simulated and real data, show that incorporating the interpretability into the model selection procedure provides more reproducible, more neurophysiologically plausible, and (as a result) more interpretable maps. Furthermore, in comparison with a standard univariate analysis, we show the proposed paradigm offers more sensitivity while preserving the interpretability of results.

Notation and Background
Let X ∈ R p be a manifold in Euclidean space that represents the input space and Y ∈ R be the output space, where Y = * (X ). Then, let S = {Z = (X, Y) | z 1 = (x 1 , y 1 ), . . . , z n = (x n , y n )} be a training set of n independently and identically distributed (i.i.d) samples drawn from the joint distribution of Z = X ×Y based on an unknown Borel probability measure ρ. In the neuroimaging context, X indicates the trials of brain recording, e.g., fMRI, MEG, or EEG signals, Y represents the experimental conditions or dependent variables, and we have S : X → Y (note the difference between S and * ). The goal of brain decoding is to find the functionˆ : X → Y as an estimation of S . Here on we refer toˆ as a brain decoding model.
As is a common assumption in the neuroimaging context, we assume the true solution of a brain decoding problem is among the family of linear functions H ( * ∈ H). Therefore, the aim of brain decoding reduces to finding an empirical approximation of S , indicated byˆ , among all ∈ H. This approximation can be obtained by estimating the predictive conditional density ρ(Y | X) by training a parametric model ρ(Y | X, ) (i.e., a likelihood function), where denotes the parameters of the model. Alternatively, can be estimated by solving a risk minimization problem: whereˆ is the parameter ofˆ , L : Z × Z → R + is the loss function, : R p → R + is the regularization term, and λ is a hyper-parameter that controls the amount of regularization. There are various choices for , each of which reduces the hypothesis space H to H ′ ⊂ H by enforcing different prior functional or structural constraints on the parameters of the linear decoding model (see for example, Tibshirani, 1996b;Tibshirani et al., 2005;Zou and Hastie, 2005;Jenatton et al., 2011). The amount of regularization λ is generally decided using cross-validation or other data perturbation methods in the model selection procedure.
In the neuroimaging context, the estimated parameters of a linear decoding modelˆ can be used in the form of a brain map so as to visualize the discriminative neurophysiological effect. Although the magnitude ofˆ (i.e., the 2-norm ofˆ ) is affected by the dynamic range of data and the level of regularization, it has no effect on the predictive power and the interpretability of maps. On the other hand, the direction ofˆ affects the predictive power and contains information regarding the importance of and relations among predictors. This type of relational information is very useful when interpreting brain maps in which the relation between different spatio-temporal independent variables can be used to describe how different brain regions interact over time for a certain cognitive process. Therefore, we refer to the normalized parameter vector of a linear brain decoder in the unit hypersphere as a multivariate brain map (MBM); we denote it by where = 2 ( . 2 represents the 2-norm of a vector). As shown in Equation (1), learning occurs using the sampled data. In other words, in the learning paradigm, we attempt to minimize the loss function with respect to S (and not * ) (Cucker and Smale, 2002). Therefore, all of the implicit assumptions (such as linearity) regarding * might not hold on S , and vice versa. The irreducible error ε is the direct consequence of sampling; it sets a lower bound on the error, where we have: The distribution of ε dictates the type of loss function L in Equation (1). For example, assuming a Gaussian distribution with mean 0 and variance σ 2 for ε implies the least squares loss function (Wu et al., 2006).

Interpretability of Multivariate Brain Maps: Theoretical Definition
In this section, we present a theoretical definition for the interpretability of linear brain decoding models and their associated MBMs. Consider a linearly separable brain decoding problem in an ideal scenario where ε = 0 and rank(X) = p.
In this case, the ideal solution of brain decoding, * , is linear and its parameters * are unique and neurophysiologically plausible (van Ede and Maris, 2016). The unique parameter vector * can be computed as follows: where X represents the covariance of X. Using * as the reference, we define the strong-interpretability of an MBM as follows: Definition 1. An MBM ˆ associated with a linear brain decoding modelˆ is "strongly-interpretable" if and only if ˆ ∝ * .
It can be shown that, in practice, the estimated solution of a linear brain decoding problem is not strongly-interpretable because of the inherent limitations of neuroimaging data, such as uncertainty (Aggarwal and Yu, 2009) in the input and output space (ε = 0), the high dimensionality of data (n ≪ p), and the high correlation between predictors (rank(X) < p). With these limitations in mind, even though in practice the solution of linear brain decoding is not strongly-interpretable, one can argue that some are more interpretable than others. For example, in the case in which * ∝ [0, 1] T , a linear classifier where ˆ ∝ [0.1, 1.2] T can be considered more interpretable than a linear classifier where ˆ ∝ [2, 1] T . This issue raises the following question: Problem 1. Let S be a training set of n i.i.d samples drawn from the joint distribution of Z = X × Y, and P(S) be the probability of drawing a certain S from Z. Assume ˆ is the MBM of a linear brain decoding modelˆ on S (estimated using Equation 1 for a certain loss function L, regularization term , and hyper-parameter λ). How can we quantify the proximity ofˆ to the strongly-intrepretable solution of the brain decoding problem * ?
To answer this question, considering the uniqueness and the plausibility of * as the two main characteristics that convey its strong-interpretability, we define the interpretability as follows: Definition 2. Let S, P(S), and ˆ be as defined in Problem 1. Then, assume α be the angle between ˆ and * . The "interpretability" (0 ≤ η ≤ 1) of a linear brain decoding modelˆ is defined as follows: In practice, only a limited number of samples are available. Therefore, perturbation techniques are used to imitate the sampling procedure. Let S 1 , . . . , S m be m perturbed training sets drawn from S via a certain perturbation scheme such as jackknife, bootstrapping (Efron, 1992), or cross-validation (Kohavi, 1995). Assume ˆ 1 , . . . , ˆ m are m MBMs estimated on the corresponding perturbed training sets, and α j (j = 1, . . . , m) be the angle between ˆ j and * . Then, the empirical version of Equation (4) can be rewritten as follows: Empirically, the interpretability is the mean of cosine similarities between * and MBMs derived fro different samplings of the training set (see Figure 1A for a schematic illustration). In addition to the fact that employing cosine similarity is a common method for measuring the similarity between vectors, we have another strong motivation for this choice. It can be shown that, for large values of p, the distribution of the dot product in the unit hyper-sphere, i.e., the cosine similarity, converges to a normal distribution with 0 mean and variance of 1 p , i.e., N (0, 1/p). Due to the small variance for a large enough p values, any similarity value that is significantly larger than zero represents a meaningful similarity between two high dimensional vectors (see Appendix 6.3 for the mathematical demonstration).
In what follows, we demonstrate how the definition of interpretability is geometrically related to the uniqueness and plausibility characteristics of the true solution of the brain decoding problem.

Interpretability Decomposition into Reproducibility and Representativeness
The trustworthy and informativeness of decoding models are providing two important motivations for improving the interpretability of models (Lipton et al., 2016). The trust of a learning algorithm refers to its ability to converge to a unique solution. On the other hand, the informativeness refers to the level of plausible information that can be derived from a model to assist or advise to a human expert. Therefore, it is expected that the interpretability can be quantified alternatively by assessing its uniqueness and neurophysiological plausibility. In this section, we firstly define the reproducibility and representativeness as measures for quantifying the uniqueness and neurophysiological plausibility of a brain decoding model, respectively. Then we show how these definitions are related to the definition of interpretability.
The high dimensionality and the high correlations between variables are two inherent characteristics of neuroimaging data that negatively affect the uniqueness of the solution of a brain decoding problem. Therefore, a certain configuration of hyper-parameters may result different estimated parameters on different portions of data. Here, we are interested in assessing this variability as a measure for uniqueness. We first define the main multivariate brain map as follows: Definition 3. Let S, P(S), and ˆ be as defined in Problem 1. The "main multivariate brain map" µ ∈ R p of a linear brain decoding modelˆ is defined as: Assuming θ j i be the ith (i = 1, . . . , p) element of an MBM estimated on the jth (j = 1, . . . , m) perturbed training set, µ empirically can be estimated by summing up ˆ j s (computed on the perturbed training set S j ) in the unit hyper-sphere: µ provides a reference for quantifying the reproducibility of an MBM: Definition 4. Let S, P(S), and ˆ be as defined in Problem 1, and µ be the main multivariate brain map ofˆ . Then, assume α be the angle between ˆ j and µ . The "reproducibility" ψ (0 ≤ ψ ≤ 1) of a linear brain decoding modelˆ is defined as follows: Let ˆ 1 , . . . , ˆ m are m MBMs estimated on the corresponding perturbed training sets, and α j (j = 1, . . . , m) be the angle between ˆ j and µ . Then, the empirical version of Eq. 8 can be rewritten as follows: In fact, reproducibility provides a measure for quantifying the dispersion of MBMs, computed over different perturbed training sets, from the main multivariate brain map. Figure 1B shows a schematic illustration for the reproducibility of a linear brain decoding model. On the other hand, the similarity between the main multivariate brain map of a decoder and the true solution can be employed as a measure for the neurophysiological plausibility of a model. We refer to this similarity as the representativeness of a linear brain decoding model: Definition 5. Let µ be the main multivariate brain map ofˆ . The "representativeness" β (0 ≤ β ≤ 1) of a linear brain decoding modelˆ is defined as the cosine similarity between its main multivariate brain map ( µ ) and the parameters of the true solution ( * ) (see Figure 1C): As discussed before, the notion of interpretabilty is tightly related to the uniqueness and plausibility, thus to the reproducibility and representativeness, of a decoding model. The following proposition analytically shows this relationship: See Appendix 6.1 for a proof. Proposition 1 indicates the interpretability of a linear brain decoding model can be decomposed into its representativeness and reproducibility. Figure 1D illustrates how the reproducibility and the representativeness of a decoding model independently affect its interpretability. Each colored region schematically represents a span of different solutions of the a certain linear model (for example with a certain configuration for its hyper-parameters) on different perturbed training sets. The area of each region schematically visualizes the reproducibility of each model, i.e., the less is the area, the higher is the reproducibility of a model. Further, the angular distance between the centroid of each region ( µ ) and the true solution ( * ) visualizes the representativeness of each corresponding model. While 1 and

A Heuristic for Practical Quantification of Interpretability in Time-Locked Analysis of MEG Data
In practice, it is impossible to evaluate the interpretability, as the true solution of the brain decoding problem * is unknown. In this study, to provide a practical proof of theoretical concepts, we exemplify contrast event-related field (cERF) as a neurophysiological plausible heuristic for the true parameters of the linear brain decoding problem ( * ) in a binary MEG decoding scenario in time domain. Due to the nature of proposed heuristic, its application is limited to the brain responses that are time-locked to the stimulus onset, i.e., the evoked responses.
The MEEG data are a mixture of several simultaneous stimulus-related and stimulus-unrelated brain activities. Assessing the electro/magneto-physiological changes that are time-locked to events of interest is a common approach to the study of MEEG data. In general, unrelated-stimulus brain activities are considered as Gaussian noise with zero mean and variance σ 2 . One popular approach to canceling the noise component is to compute the average of multiple trials. The assumption is that, when the effect of interest is time-locked to the stimulus onset, the independent noise components can be vanished by means of averaging. It is expected that the average will converge to the true value of the signal with a variance of σ 2 n (where n is the number of trials). The result of the averaging process consist of a series of positive and negative peaks occurring at a fixed time relative to the event onset, generally known as ERF in the MEG context. These component peaks are reflecting phasic activity that are indexed with different aspects of cognitive processing (Rugg and Coles, 1995) 1 .
Assume X + = {x i ∈ X | y i = 1} ∈ R n + ×p and X − = {x i ∈ X | y i = −1} ∈ R n − ×p be sets of positive and negative samples in a binary MEG decoding scenario. Then, the cERF brain map cERF is computed as follows: Generally speaking cERF is a contrast ERF map between two experimental conditions. Using the core theory presented in Haufe et al. (2013), the equivalent generative model for the solution of linear brain decoding, i.e., the activation pattern (A), is unique and we have: Assumingˆ to be the solution of least squares in a binary decoding scenario, then the following proposition describes the relation between cERF and the activation pattern A: See Appendix 6.2 for the proof. Proposition 2.4 shows that, in a binary time-domain MEG decoding scenario, cERF is proportional to the equivalent generative model for the solution of least squares classifier. Furthermore, cERF is proportional to the t-statistic that is widely used in the univariate analysis of neuroimaging data. Using cERF as a heuristic for * , the representativeness can be approximated as follows: Whereβ is an approximation of the actual representativeness β . In a similar manner, cERF can be used to heuristically approximate the interpretability as follows: where γ 1 , . . . , γ m are the angles between ˆ 1 , . . . , ˆ m and cERF . It can be shown thatη =β × ψ .
The proposed heuristic is only applicable to the evoked responses in sensor and source space MEEG data. Despite this limitation, cERF provides an empirical example that shows how the presented theoretical definitions can be applied in a real decoding scenario. The choice of the heuristic has a direct effect on the approximation of interpretability and that an inappropriate selection of the heuristic yields a very poor estimation of interpretability. Therefore, the choice of heuristic should be carefully justified based on accepted and well-defined facts regarding the nature of the collected data.
Since the labels are used in the computation of cERF, a proper validation strategy should be employed to avoid the double dipping issue (Kriegeskorte et al., 2009). One possible approach is to exclude the entire test set from the model selection procedure using a nested nested cross-validation strategy. An alternative approach is employing model averaging techniques to neatly get advantage of the whole dataset (Varoquaux et al., 2017). Since our focus is on the model selection, in the remaining text, we implicitly assume the test data is excluded from the experiments, thus, all the experimental results are reported on the training and validation sets.

Incorporating the Interpretability into Model Selection
The procedure for evaluating the performance of a model so as to choose the best values for hyper-parameters is known as model selection (Hastie et al., 2009). This procedure generally involves numerical optimization of the model selection criterion on the training and validation sets (and not the test set). Let U be a set of hyper-parameters, then the goal of model selection procedure reduces to finding the best model configuration u * ∈ U that maximizes the model selection criterion (e.g., generalization performance) on the training set S. The most common model selection criterion is based on an estimator of generalization performance, i.e., the predictive power. In the context of brain decoding, especially when the interpretability of brain maps matters, employing the predictive power as the only decisive criterion in model selection is problematic in terms of interpretability of MBMs Rasmussen et al., 2012;Conroy et al., 2013;Varoquaux et al., 2017). Valverde-Albacete and Peláez-Moreno (2014) experimentally showed that in a classification task optimizing only classification error rate is insufficient to capture the transfer of crucial information from the input to the output of a classifier. This fact highlights the importance of having some control over the estimated model weights in the model selection. Here, we propose a multiobjective criterion for model selection that takes into account both prediction accuracy and MBM interpretability.
Letη and δ be the approximated interpretability and the generalization performance of a linear brain decoding modelˆ , respectively. We propose the use of the scalarization technique (Caramia and Dell´Olmo, 2008) for combiningη and δ into one scalar 0 ≤ ζ ( ) ≤ 1 as follows: where ω 1 and ω 2 are weights that specify the level of importance of the interpretability and the performance, respectively. κ is a threshold on the performance that filters out solutions with poor performance. In classification scenarios, κ can be set by adding a small safe interval to the chance level of classification.
The hyper-parameters that are optimized based on ζ are Pareto optimal (Marler and Arora, 2004). We hypothesize that optimizing the hyper-parameters based on ζ , rather only δ , yields more informative MBMs. Algorithm 1 summarizes the proposed model selection scheme. The model selection procedure receives the training set S and a set of possible configurations for hyper-parameters U, and returns the best hyper-parameter configuration u * .

Toy Dataset
We regenerate the simple 2-dimensional toy data presented in Haufe et al. (2013). Assume that the true underlying generative function * is defined by: [−1.5, 0] T }; and x 1 and x 2 represent the first and the second dimension of the data, respectively. Furthermore, assume the data is contaminated by Gaussian noise with co-variance = 1.02 −0.3 −0.3 0.15 . In fact, the Gaussian noise adds uncertainty to the input space.

Simulated MEG Data
We simulated two classes of MEG data, each of which composed of 250 epochs with length of 330ms at 300Hz sampling rate (so that we have 100 time-points). For simplicity, the whole scalp topography are simulated with a single dipole located at −4.7, −3.7, and 5.3cm in the RAS (right, anterior, superior) coordinate system. The dipole is oriented toward [1,1,0] direction in the RA plane (see Figure 2A). One hundred two magnetometer sensors of Elekta Neuromag system are simulated using a standard forward model algorithm implemented in the Fieldtrip toolbox (Oostenveld et al., 2010). Compute cERF on S. ⊲ using Equation (11) 3: for all u i ∈ U do ⊲ For all hyper-parameter configurations.

5:
Partition S into training S tr and validation S vl 6: subsets via a perturbation method.

7:
Computeˆ j on S tr using u i as the 8: hyper-parameter. end 9: Compute δ i ofˆ j s on S vl .
⊲ using Equation (14) 11: Compute ζ i . ⊲ using Equation (15) end 12: u * = argmax u i ∈U (ζ ). 13: return u * . time-course: (1) a time-locked ERF effect with a positive 3Hz followed by a negative 5 Hz half-cycle sinusoid peaks after 150 ± 10ms and 250 ± 10ms of the epoch onset, respectively; (2) uncorrelated background brain activity that was simulated by summing 50 sinusoids with random frequency from 1 to 125Hz, and random phase varied between 0 and 2π. Following the data simulation procedure in Yeung et al. (2004), the amplitude of any single frequency component of the signal (the ERF effect and the background noise) is set based on the empirical spectral power of human brain activity to mimic the actual magnetic features of scalp surface; and (3) white Gaussian noise scaled with the root mean squared of the signal in each epoch. The epochs of the negative class are constructed without the ERF effect by adding up only the noise components (i.e., the background activity and the white noise). Therefore, the ERF component is considered as the discriminative ground-truth in our experiments (see Figure 2B).

MEG Data
We use the MEG dataset presented in Henson et al. (2011) 2 . The dataset was also used for the DecMeg2014 competition 3 . In this dataset, visual stimuli consisting of famous faces, unfamiliar faces, and scrambled faces are presented to 16 subjects and fMRI, EEG, and MEG signals are recorded. Here, we are only interested in MEG recordings. The MEG data were recorded using a VectorView system (Elekta Neuromag, Helsinki, Finland) with a magnetometer and two orthogonal planar gradiometers located at 102 positions in a hemispherical array in a light Elekta-Neuromag magnetically shielded room.
Three major reasons motivated the choice of this dataset: (1) It is publicly available.
(2) The spatio-temporal dynamic of the MEG signal for face vs. scramble stimuli has been well studied. The event-related potential analysis of EEG/MEG shows that N170 occurs 130 − 200ms after stimulus presentation and reflects the neural processing of faces (Bentin et al., 1996;Henson et al., 2011). Therefore, the N170 component can be considered the ground truth for our analysis. (3) In the literature, non-parametric mass-univariate analysis such as cluster-based permutation tests is unable to identify narrowly distributed effects in space and time (e.g., an N170 component; Groppe et al., 2011a,b). These facts motivate us to employ multivariate approaches that are more sensitive to these effects.
As in Olivetti et al. (2014), we created a balanced face vs. scrambled MEG dataset by randomly drawing from the trials of unscrambled (famous or unfamiliar) faces and scrambled faces in equal number. The samples in the face and scrambled face categories are labeled as 1 and −1, respectively. The raw data is high-pass filtered at 1Hz, down-sampled to 250Hz, and trimmed from 200ms before the stimulus onset to 800ms after the stimulus. Thus, each trial has 250 time-points for each of the 306 MEG

Classification and Evaluation
In all experiments, Lasso (Tibshirani, 1996b) classifier with ℓ 1 penalization is used for decoding. Lasso is a very popular classification method in the context of brain decoding, mainly because of its sparsity assumption. The choice of Lasso, as a simple model with only one hyper-parameter, helps us to better illustrate the importance of including the interpretability in the model selection (see the supplementary materials for the results of the elastic-net; Zou and Hastie, 2005 classifier). The solution of decoding is computed by solving the following optimization problem:ˆ = argmin L(X , Y) + λ 1 (16) where . 1 represents the ℓ 1 -norm, and λ is the hyper-parameter that specifies the level of regularization. Therefore, the aim of the model selection is to find the best value for λ on the training set S. Here, we try to find the best regularization parameter value among λ = {0.001, 0.01, 0.1, 1, 10, 50, 100, 250, 500, 1000}.
We use the out-of-bag (OOB) (Wolpert and Macready, 1999;Breiman, 2001) method for computing δ , ψ ,β ,η , and ζ for different values of λ. In OOB, given a training set (X, Y), m replications of bootstrap (Efron, 1992) are used to create perturbed training and validation sets (we set m = 50) 5 . In all of our experiments, we set ω 1 = ω 2 = 1 and κ = 0.6 in the computation of ζ . Furthermore, we set δ = 1−EPE where EPE indicates the expected prediction error; it is computed using the procedure explained in Appendix 6.4. Employing OOB provides the possibility of computing the bias and variance of the model as contributing factors in EPE.

Performance-Interpretability Dilemma: A Toy Example
In the definition of * on the toy dataset discussed in Section 2.6.1, x 1 is the decisive variable and x 2 has no effect on the classification of samples into target classes. Therefore, excluding the effect of noise and based on the theory of the maximal margin classifier (Vapnik and Kotz, 1982), * ∝ [1, 0] T is the true solution to the decoding problem. By accounting for the effect of noise, solving the decoding problem in (X, Y) space yields ˆ ∝ [1/ √ 5, 2/ √ 5] T as the parameters of the linear classifier. Although the estimated parameters on the noisy data provide the best generalization performance for the noisy samples, any attempt to interpret this solution fails, as it yields the wrong conclusion with respect to the ground truth (it says x 2 has twice the influence of x 1 on the results, whereas it has no effect). This simple experiment shows that the most accurate model is not always the most interpretable one, primarily because the contribution of the noise in the decoding process (Haufe et al., 2013). On the other hand, the true solution of the problem * does not provide the best generalization performance for the noisy data.
To illustrate the effect of incorporating the interpretability in the model selection, a Lasso model with different λ values is used for classifying the toy data. In this example, because * is known, the exact value of interpretability can be computed using Equation (5). Table 1 compares the resultant performance and interpretability from Lasso. Lasso achieves its highest performance (δ = 0.9884) at λ = 10 with ˆ ∝ [0.4636, 0.8660] T (indicated by the black dashed line in Figure 3). Despite having the highest performance, this solution suffers from a lack of interpretability (η = 0.4484). By increasing λ, the interpretability improves so that for λ = 500, 1000 the classifier reaches its highest interpretability by compensating for 0.06 of its performance. Our observation highlights two main points: 1. In the case of noisy data, the interpretability of a decoding model can be possibly incoherent with its performance. Thus, optimizing the parameter of the model based on its performance does not necessarily improve its interpretability. This observation confirms the previous finding by Rasmussen et al. (2012) regarding the tradeoff between the spatial reproducibility (as a measure for the interpretability) and the prediction accuracy in brain decoding. 2. If the right criterion is used in the model selection, employing proper regularization technique (sparsity prior, in the case of toy data) leads to more interpretable decoding models.

Decoding on Simulated MEG Data
With the main aim of comparing the quality of the heuristically approximated interpretability with respect to its actual value, we solve the decoding problem on the simulated MEG data where the ground-truth discriminative effect is known. The ground truth effect * is used to compute the actual interpretability of the decoding model. On the other hand, interpretability is approximated by means of cERF . The whole data simulation and decoding processes are repeated 25 times and the results are summarized in Figure 4. Figures 4A,B show the actual   (η ) and the approximated (η ) interpretability for different λ values. Even thoughη consistently overestimates η , there is a significant co-variation (Pearson's correlation p-value = 9 × 10 −4 ) between two measures as λ increases. Thus, despite overestimation problem of the heuristically approximated interpretability values, they are still reliable measures for quantitative comparison between interpretability level of brain decoding models with different hyper-parameters. For example, both η andη suggest the decoding model with λ = 50 as the most interpretable model. Figure 4C shows brain decoding models at λ = 10 and λ = 50 yield equivalent generalization performances (Wilcoxon rank sum test p-value = 0.08), while the MBM resulted from λ = 50 has significantly higher interpretability (Wilcoxon rank sum test p-value = 4 × 10 −9 ). The advantage of this difference in interpretability levels is visualized in Figure 5 where topographic maps are plotted for the weights of brain decoding models with different λ values by averaging the classifier weights in the time interval of 100-200 ms. The visual comparison shows MBM at λ = 50 is more similar to the ground-truth map (see Figure 2B) than the MBMs computed at other λ values. This superiority is well-reflected in the corresponding approximated interpretability values, that confirms the effectiveness of the interpretability criterion in measuring the level of information in the MBMs.
The results of this experiment confirm again the fact that the generalization performance is not a reliable criterion to measure the level of information learned by a linear classifier. For example consider the decoding model with λ = 1 in which the performance of the model is significantly above the chance level (see Figure 4C) while the corresponding MBM ( Figure 5A) is completely misrepresents the ground-truth effect (Figure 2B).

Single-Subject Decoding on MEG Data
To investigate the behavior of the proposed model selection criterion ζ , we benchmark it against the commonly used performance criterion δ in a single-subject decoding scenario. Assuming (X i , Y i ) for i = 1, . . . , 16 are MEG trial/label pairs for subject i, we separately train a Lasso model for each subject to estimate the parameter of the linear functionˆ i , where Y i = X iˆ i . We represent the optimized solution based on δ and ζ byˆ δ i andˆ ζ i , respectively. We also denote the MBM associated withˆ δ i andˆ ζ i by ˆ δ i and ˆ ζ i , respectively. Therefore, for each subject, we compare the resulting decoders and MBMs computed based on these two model selection criteria.  Frontiers in Neuroscience | www.frontiersin.org Figure 6A represents the mean and standard-deviation of the performance and interpretability of Lasso across 16 subjects for different λ values. The performance and interpretability curves further illustrate the performance-interpretability dilemma of Lasso classifier in the single-subject decoding scenario, in which increasing the performance delivers less interpretability. The average performance across subjects is improved when λ approaches 1, but on the other side, the reproducibility and the representativeness of models declines significantly (see Figure 6B; Wilcoxon rank sum test p-value = 9 × 10 −4 and 8 × 10 −7 , respectively). In fact, in this dataset a higher amount of sparsity increases the gap between the generalization performance and interpretability.
One possible reason behind the performance-interpretability dilemma in this experiment is illustrated in Figure 6C. The figure shows the mean and standard deviation of bias, variance, and EPE of Lasso across 16 subjects. The plot shows while the change in bias is correlated with that of EPE (Pearson's correlation coefficient = 0.9993), there is anti-correlation between the trends of variance and EPE (Pearson's correlation coefficient = −0.8884). Furthermore, it proposes that the effect of variance is overwhelmed by bias in the computation of EPE, where the best performance (minimum EPE) at λ = 1 has the lowest bias, its variance is higher than for λ = 0.001, 0.01, 0.1. While this tiny increase in the variance has negligible effect on the EPE of the model, Figure 6B shows its significant (Wilcoxon rank sum test p-value = 8 × 10 −7 ) negative effect on the reproducibility of maps from λ = 0.1 to λ = 1. Table 2 summarizes the performance, reproducibility, representativeness, and interpretability ofˆ δ i andˆ ζ i for 16 subjects. The average result over 16 subjects shows that employing ζ instead of δ in model selection provides higher reproducibility, representativeness, and (as a result) interpretability compensating for 0.04 of performance. The last column of table (δ cERF ) summarizes the performance of decoding models over 16 subjects when cERF is used as classifier weights. The comparison illustrates a significant difference (Wilcoxon rank sum test p-value = 1.5 × 10 −6 ) between δ cERF and δ( )s. These facts demonstrate that ˆ ζ is a good compromise between ˆ δ and cERF in terms of classification performance and model interpretability.
These results are further analyzed in Figure 7 whereˆ δ i and ζ i are compared subject-wise in terms of their performance and interpretability. The comparison shows that adopting ζ instead of δ as the criterion for model selection yields higher interpretable models by compensating a negligible degree of performance in 14 out of 16 subjects. Figure 7A shows that employing δ provides on average slightly higher accurate models (Wilcoxon rank sum test p-value = 0.012) across subjects (0.83 ± 0.05) than using ζ (0.79 ± 0.04). On the other side, Figure 7B shows that employing ζ and compensating by 0.04 in the performance provides (on average) substantially higher (Wilcoxon rank sum test p-value= 5.6 × 10 −6 ) interpretability across subjects (0.62 ± 0.05) compared to δ (0.31 ± 0.12). For example, in the case of subject 1 (see Table 2), using δ in model selection to select the best λ value for the Lasso yields a model with δ = 0.81 andη = 0.26. In contrast, using ζ delivers a model with δ = 0.78 andη = 0.63. This inverse relationship between performance and interpretability is direct consequence of over-fitting of model to the noise structure in a small-sample size brain decoding problem (see Section 3.1).
The advantage of the exchange between the performance and the interpretability can be seen in the quality of MBMs. Figures 8A,B show ˆ δ 1 and ˆ ζ 1 of subject 1, i.e., the spatiotemporal multivariate maps of the Lasso models with maximum values of δ and ζ , respectively. The maps are plotted for 102 magnetometer sensors. In each case, the time course of weights of classifiers associated with the MEG2041 and MEG1931 sensors are plotted. Furthermore, the topographic maps represent the spatial patterns of weights averaged between 184ms and 236ms after stimulus onset. While ˆ δ 1 is sparse in time and space, it fails to accurately represent the spatio-temporal dynamic of the N170 component. Furthermore, the multicollinearity problem arising from the correlation between the time course of the MEG2041 and MEG1931 sensors causes extra attenuation of the N170 effect   The bold indicates the best mean values over different criteria. in the MEG1931 sensor. Therefore, the model is unable to capture the spatial pattern of the dipole in the posterior area. In contrast, ˆ ζ 1 represents the dynamic of the N170 component in time.
In addition, it also shows the spatial pattern of two dipoles in the posterior and temporal areas. In summary, ˆ ζ 1 suggests a more representative pattern of the underlying neurophysiological effect than ˆ δ 1 . In addition, optimizing the hyper-parameters of brain decoding based on ζ offers more reproducible brain decoders. According to Table 2, using ζ instead of δ provides (on average) 0.15 more reproducibility over 16 subjects. To illustrate the advantage of higher reproducibility on the interpretability of maps, Figure 9 visualizes ˆ δ 1 and ˆ ζ 1 over 4 perturbed training sets. The spatial maps (Figures 9A,C) are plotted for the magnetometer sensors averaged in the time interval between 184ms and 236ms after stimulus onset. The temporal maps (Figures 9B,D) are showing the multivariate temporal maps of MEG1931 and MEG2041 sensors. While ˆ δ 1 is unstable in time and space across the 4 perturbed training sets, ˆ ζ 1 provides more reproducible maps.

Mass-Univariate Hypothesis Testing on MEG Data
It is shown by Groppe et al. (2011a,b) that non-parametric mass-univariate analysis is unable to detect narrowly distributed effects in space and time (e.g., an N170 component). To 1 provides a better spatio-temporal representation of the N170 effect than ˆ δ 1 .
illustrate the advantage of the proposed decoding framework for spotting these effects, we performed a non-parametric cluster-based permutation test (Maris and Oostenveld, 2007) on our MEG dataset using Fieldtrip toolbox (Oostenveld et al., 2010). In a single subject analysis scenario, we considered the trials of MEG recordings as the unit of observation in a between-trials experiment. Independent-samples t-statistics are used as the statistics for evaluating the effect at the sample level and to construct spatio-temporal clusters. The maximum of the cluster-level summed t-value is used for the cluster level statistics; the significance probability is computed using a Monte Carlo method. The minimum number of neighboring channels for computing the clusters is set to 2. Considering 0.025 as the two-sided threshold for testing the significance level and repeating the procedure separately for magnetometers and combined-gradiometers, no significant result is found for any of the 16 subjects. This result motivates the search for more sensitive (and, at the same time, more interpretable) alternatives for univariate hypothesis testing.

Defining Interpretability: Theoretical Advantages
An overview of the brain decoding literature shows frequent co-occurrence of the terms interpretation, interpretable, and interpretability with the terms model, classification, parameter, decoding, method, feature, and pattern (see the quick metaanalysis on the literature in the supplementary material); however, a formal formulation of the interpretability is never presented. In this study, our primary interest is to present a simple and theoretical definition of the interpretability of linear brain decoding models and their corresponding MBMs. Furthermore, we show the way in which interpretability is related to the reproducibility and neurophysiological representativeness of MBMs. Our definition and quantification of interpretability remains theoretical, as we assume that the true solution of the brain decoding problem is available. Despite this limitation, we argue that the presented definition provides a concrete framework of a previously abstract concept and that it establishes a theoretical background to explain an ambiguous phenomenon in the brain decoding context. We support our argument using an example in the time-domain MEG decoding in which we show how the presented definition can be exploited to heuristically approximate the interpretability. Our experimental results on MEG data shows accounting for the approximated measure of interpretability has a positive effect on the human interpretation of brain decoding models. This example shows how partial prior knowledge regarding the timing and location of neural activity can be used to find more plausible multivariate patterns in data. Furthermore, the proposed decomposition of the interpretability of MBMs into their reproducibility and representativeness explains the relationship between the influential cooperative factors in the interpretability of brain decoding models and highlights the possibility of indirect and partial evaluation of interpretability by measuring these effective factors.

Application in Model Evaluation
Discriminative models in the framework of brain decoding provide higher sensitivity and specificity than univariate analysis in hypothesis testing of neuroimaging data. Although multivariate hypothesis testing is performed based solely on the generalization performance of classifiers, the emergent need for extracting reliable complementary information regarding the underlying neuronal activity motivated a considerable amount of research on improving and assessing the interpretability of classifiers and their associated MBMs. Despite ubiquitous use, the generalization performance of classifiers is not a reliable criterion for assessing the interpretability of brain decoding models (Rasmussen et al., 2012;Varoquaux et al., 2017). Therefore, considering extra criteria might be required. However, because of the lack of a formal definition for interpretability, different characteristics of linear classifiers are considered as the decisive criterion in assessing their interpretability. Reproducibility (Rasmussen et al., 2012;Conroy et al., 2013), stability selection Wang et al., 2015), sparsity (Dash et al., 2015;Shervashidze and Bach, 2015), and neurophysiological plausibility (Afshin-Pour et al., 2011) are examples of related criteria.
Our definition of interpretability helped us to fill this gap by introducing a new multi-objective model selection criterion as a weighted compromise between interpretability and generalization performance of linear models. Our experimental results on single-subject decoding showed that adopting the new criterion for optimizing the hyper-parameters of brain decoding models is an important step toward reliable visualization of learned models from neuroimaging data. It is not the first time in the neuroimaging context that a new metric is proposed in combination with generalization performance for the model selection. Several recent studies proposed the combination of the reproducibility of the maps (Rasmussen et al., 2012;Conroy et al., 2013;Strother et al., 2014) or the stability of the classifiers (Yu, 2013;Lim and Yu, 2016;Varoquaux et al., 2017) with the performance of discriminative models to enhance the interpretability of decoding models. Our definition of interpretability supports the claim that the reproducibility is not the only effective factor in interpretability. Therefore, our contribution can be considered a complementary effort with respect to the state of the art of improving the interpretability of brain decoding at the model selection level.
Furthermore, this work presents an effective approach for evaluating the quality of different regularization strategies for improving the interpretability of MBMs. As briefly reviewed in Section 1, there is a trend of research within the brain decoding context in which the prior knowledge is injected into the decoding process via the penalization term in order to improve the interpretability of decoding models. Thus far, in the literature, there is no ad-hoc method to directly compare the interpretability of MBMs resulting from different penalization techniques. Our findings provide a further step toward direct evaluation of interpretability of the currently proposed penalization strategies. Such an evaluation can highlight the advantages and disadvantages of applying different strategies on different data types and facilitates the choice of appropriate methods for a certain application. Haufe et al. (2013) demonstrated that the weight in linear discriminative models are unable to accurately assess the relationship between independent variables, primarily because of the contribution of noise in the decoding process. The authors concluded that the interpretability of brain decoding cannot be improved using regularization. The problem is primarily caused by the decoding process per se, where it minimizes the classification error only considering the uncertainty in the output space (Zhang, 2005;Aggarwal and Yu, 2009;Tzelepis et al., 2015) and not the uncertainty in the input space (or noise). Our experimental results on the toy data (see Section 3.1) shows that if the right criterion is used for selecting the best values for hyperparameters, appropriate choice of the regularization strategy can still play a significant role in improving the interpretability of results. For example, in the case of toy data, the true generative function behind the sampled data is sparse (see Section 2.6.1), but because of the noise in the data, the sparse model is not the most accurate one. On the other hand, a more comprehensive criterion (in this case, ζ ) that considers also the interpretability of model parameters facilitates the selection of correct prior assumptions about the distribution of the data via regularization. This observation encourages the modification of the conclusion in Haufe et al. (2013) as follows: if the performance of the model is the only criterion in the model selection, then the interpretability cannot necessarily be improved by means of regularization. This modification offers a practical shift in methodology, where we propose to replace the post-processing of weights proposed in Haufe et al. (2013) with refinement of hyper-parameter selection based on the newly developed model selection criterion.

The Performance-Interpretability Dilemma
The performance-interpretability dilemma refers to the trade-off between the generalization performance and the interpretability of a decoding model. In some applications of brain decoding, such as BCI, a more accurate model (even with no interpretability) is desired. On the other hand, when the brain decoding is employed for hypothesis testing purpose, an astute balance between two factors is more favorable. The presented metric for model selection (ζ ) provides the possibility to maintain this balance. An important question at this point is on the nature of the performance-interpretability dilemma, whether it is model-driven or data-driven? In other words, whether some decoding models (e.g., sparse models) suffer from this deficit, or it is independent from the decoding model and depends on the distribution of data rather assumptions of the decoding model.
Our experiments shed light on the fact that the performanceinterpretability dilemma is driven by the uncertainty (Aggarwal and Yu, 2009) in data. The uncertainty in data refers to the difference between the true solution of decoding * and the solution of decoding in sampled data space S , and is generally consequence of noise in the input or/and output spaces. This gap between * and S is also known as irreducible error (see Equation 2) in the learning theory, and it cannot fundamentally be reduced by minimizing the error. Therefore, any attempt toward improving the classification performance in the sampled data space might increase the irreducible error. As an example, our experiment on the toy data (see Section 3.1) shows the effect of noise in input space on the performance-interpretability dilemma. Improving the performance of the model (i.e., fitting to S ) diverges the estimated solution of decodingˆ from its true solution * , thus reduces the interpretability of the decoding model. Furthermore, our experiments demonstrate that incorporating the interpretability of decoding models in model selection facilitates finding the best match between the decoding model and the distribution of data. For example in classification of toy data, the new model selection metric ζ selects the more sparse model with a better match to the true distribution of data, despite worse generalization performance.

Advantage over Mass-Univariate Analysis
Mass-univariate hypothesis testing methods are among the most popular tools for forward inference on neuroimaging data in cognitive neuroscience field. Mass-univariate analyses consist of univariate statistical tests on single independent variables followed by multiple comparison correction. Generally, multiple comparison correction reduces the sensitivity of massunivariate approaches because of the large number of univariate tests involved. Cluster-based permutation testing (Maris and Oostenveld, 2007) provides a more sensitive univariate analysis framework by making the cluster assumption in the multiple comparison correction. Unfortunately, this method is not able to detect narrow spatio-temporal effects in the data (Groppe et al., 2011a). As a remedy, brain decoding provides a very sensitive tool for hypothesis testing; it has the ability to detect multivariate patterns, but suffers from a low level of interpretability. Our study proposes a possible solution for the interpretability problem of classifiers, and therefore, it facilitates the application of brain decoding in the analysis of neuroimaging data. Our experimental results for the MEG data demonstrate that, although the nonparametric cluster-based permutation test is unable to detect the N170 effect in MEG data, employing ζ instead of δ in model selection not only detects the stimuli-relevant information in the data, but also assures both reproducible and representative spatio-temporal mapping of the timing and the location of underlying neurophysiological effect.

Limitations and Future Directions
Despite theoretical and practical advantages, the proposed definition and quantification of interpretability suffer from some limitations. All of the presented concepts are defined for linear models, with the main assumption that * ∈ H (where H is a class of linear functions). This fact highlights the importance of linearizing the experimental protocol in the data collection phase (Naselaris et al., 2011). Extending the definition of interpretability to non-linear models demands future research into the visualization of non-linear models in the form of brain maps. Currently, our findings cannot be directly applied to non-linear models. Furthermore, the proposed heuristic for the time-domain MEG data applies only to binary classification. One possible solution in multiclass classification is to separate the decoding problem into several binary sub-problems. In addition the quality of the proposed heuristic is limited for the small sample size datasets. Of course the proposed heuristic is just an example of possible options for assessing the neurophysiological plausibility of MBMs in time-locked analysis of MEG data, thus, improving the quality of heuristic would be of interest in future researches. Finding physiologically relevant heuristics for other acquisition modalities such as fMRI, or frequency domain MEEG data, can be also considered as possible directions in future work.

CONCLUSIONS
We presented a novel theoretical definition for the interpretability of linear brain decoding and associated multivariate brain maps. We demonstrated how the interpretability relates to the representativeness and reproducibility of brain decoding. Although it is theoretical, the presented definition provides a first step toward practical solution for filling the knowledge extraction gap in linear brain decoding. As an example of this major breakthrough, and to provide a proof of concept, a heuristic approach based on the contrast event-related field is proposed for practical evaluation of the interpretability in multivariate recovery of evoked MEG responses. We experimentally showed that adding the interpretability of brain decoding models as a criterion in the model selection procedure yields significantly higher interpretable models by sacrificing a negligible amount of performance. Our methodological and experimental achievements can be considered a complementary theoretical and practical effort that contributes to researches on enhancing the interpretability of multivariate pattern analysis.

AUTHOR CONTRIBUTIONS
SK contributed in developing the theoretical and experimental contents of this study. SV and AP were involved in developing the theoretical machine learning aspects. NW was advising on the MEG experimental aspects.

(A1)
A similar procedure can be used to proveη =β × ψ by replacing * with cERF .

A.2. Proof of Proposition 2
According to Haufe et al. (2013), for a linear discriminative model with parametersˆ , the unique equivalent generative model can be computed as follows: In a binary (Y = {1, −1}) least squares classification scenario, we have: where X represents the covariance of the input matrix X, and µ + and µ − are the means of positive and negative samples, FIGURE A1 | Relation between representativeness, reproducibility, and interpretability in two dimensions.
respectively. Therefore, the equivalent generative model for the above classification problem can be derived by computing the difference between the mean of samples in two classes that is equivalent to the definition of cERF in time-domain MEG data.

A.3. The Distribution of Cosine Similarity
The aim of this section is to illustrate that the probability density function (PDF) of the cosine similarity between two randomly drawn vectors in the high dimensional space (large p) is very close to normal distribution with zero mean and small variance.
To do this, we first need to find the distribution of dot product in the uniform unit hyper-sphere. Let a and b be two uniformly drawn random vectors from a unit hyper-sphere in R p . Assuming that γ is the angle between a and b, the distribution of cosine similarity is equivalent to the dot product < a.b >. Without loss of generality, let b be along the positive x-axis in the coordinate system. Thus, the dot product < a.b > is the projection of a on the x-axis, i.e., x coordinate of a. Therefore, for a certain value of γ , the dot product is a p − 1 dimensional hyper-sphere that is orthogonal to the x-axis (the red circle in Figure A2) and the PDF of the dot product is the surface area of p dimensional hypersphere constructed by the dot products for different γ values (the dashed blue sphere in Figure A2). To compute the area of this hyper-sphere we take the sum of the surface area of the p dimensional conical frustums over small intervals dx (the gray area in Figure A2): where (1 − x 2 ) p−2 is the surface area of the base of the cone (e.g., the perimeter of the red circle in Figure A2) and dx 1−x 2 is the slope