Let X∈ℝp be a manifold in Euclidean space that represents the input space and Y∈ℝ be the output space, where Y=Φ*(X). Then, let S = {Z = (X, Y) | z₁ = (x₁, y₁), …, 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→ℝ+ is the loss function, Ω : ℝ^p → ℝ⁺ 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 hyper-sphere as a multivariate brain map (MBM); we denote it by Θ→ where Θ→=Θ||Θ||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 σ² for ε implies the least squares loss function (Wu et al., 2006).

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: (4)ηΦ=𝔼P(S)[cos(α)]

In practice, only a limited number of samples are available. Therefore, perturbation techniques are used to imitate the sampling procedure. Let S¹, …, 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 1p, 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.

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” Θ→μ∈ℝp of a linear brain decoding model Φ^ is defined as: (6)Θ→μ=𝔼P(S)[Θ^→]||𝔼P(S)[Θ^→]||2

Assuming θij 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 Θ^→js (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: (8)ψΦ=EP(S)[cos(α)]

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): (10)βΦ=|Θ→μ.Θ→*|||Θ→μ||2||Θ→*||2

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:

Proposition 1. η_Φ = β_Φ × ψ_Φ.

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 Φ₁ and Φ₂ have similar reproducibility, Φ₂ has higher interpretability than Φ₁ because it is more representative of the true solution. On the other hand, Φ₁ and Φ₃ have similar representativeness, however Φ₃ is more interpretable due to the higher level of reproducibility.

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 σ². 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 σ2n (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)¹.

Assume X+={xi∈X|yi=1}∈ℝn+×p and X-={xi∈X|yi=-1}∈ℝ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:

Proposition 2. Θ→cERF∝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 γ₁, …, γ_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.

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 (Gramfort et al., 2012; 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 multi-objective 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 ω₁ and ω₂ 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^*.

In all experiments, Lasso (Tibshirani, 1996b) classifier with ℓ₁ 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:

where ||.||₁ represents the ℓ₁-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)⁵. In all of our experiments, we set ω₁ = ω₂ = 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.

In the definition of Φ^* on the toy dataset discussed in Section 2.6.1, x₁ is the decisive variable and x₂ 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₂ has twice the influence of x₁ 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:

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 trade-off between the spatial reproducibility (as a measure for the interpretability) and the prediction accuracy in brain decoding.

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⁻⁴) 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⁻⁹). 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).

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 Yi=XiΘ^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.

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⁻⁴ and 8 × 10⁻⁷, 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⁻⁷) 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⁻⁶) 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⁻⁶) 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 spatio-temporal 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 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.

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 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.

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 meta-analysis 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.

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 (Varoquaux et al., 2012; 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 hyper-parameters, 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 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 performance-interpretability 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.

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 mass-univariate 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 non-parametric 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.

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.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.