Basic definitions

All symbols used in the following derivations are summarized in Table 2. First, we define U to be the matrix of noiseless activity profiles with K (number of experimental conditions) rows and P (number of measurement channels) columns. Each row of this matrix is an activity pattern, the response of the whole population to a single condition. Each column of this matrix is an activity profile (Fig 1A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Notation used.

For non-scalars, the second column indicates the vector / matrix size.

https://doi.org/10.1371/journal.pcbi.1005508.t002

Because we are interested in the distribution of activity profiles, but not in the activity profiles per se, we consider the columns of U to be a random variable. This is an essential step underlying our common framework, which is justified by the fact that, for the purpose of reading out information, the different measurement channels are exchangeable (see introduction). We assume that the activity profiles are repeatedly measured, with the data consisting of M independent partitions, each containing at least one activity measurement for each condition and measurement channel. In the context of fMRI, a partition will consist of a separate phase of data acquisition, e.g. a scanner run. The activity estimates of partition m are the true patterns U plus noise E^(m). The noise captures both neural trial-by-trial variability of the activity pattern in a single condition, as well as measurement noise. (1)

For the purposes of this paper, we assume that the noise is Gaussian, and independent and identically distributed (i.i.d.) across conditions and partitions (homoscedasticity). Possible dependence within each partition, however, can be easily accounted for [29, 31].

Dependence between measurement channels

The discussion below further assumes that the noise is also i.i.d. across different measurement channels (isotropicity). However, noise in fMRI, MEG, and even invasive electrophysiology exhibits strong correlations between neighboring locations in the brain. To account for these dependencies, we employ multivariate noise normalization (i.e. spatial prewhitening), which has been shown to increase the reliability of inference [32]. Across all measurement channels, we estimate the P × P variance-covariance matrix across trials, Σ_P and then regularize the estimate by shrinking it towards a diagonal matrix [33]. In the context of fMRI, we can use the residual time series from the fitting of the time-series model to estimate noise covariance [32, 34]. We then post-multiply our activity estimates by , rendering the model errors in the channels approximately uncorrelated. If multivariate noise normalization is not performed or is incomplete, inference will be suboptimal in all three methods (for details see [29]).

Second moment matrix and decodeability

In this section, we show that, under the assumption of Gaussian noise, the second moment of the activity profiles fully characterizes the decodability of any feature that describes the experimental conditions. The fact that the second moment determines what can be decoded provides a motivation, from the perspective of brain computation, for using the second moment matrix as a summary statistic.

The n^th moment of a scalar random variable u is E(uⁿ), where E() denotes the expected value. Here we use a multivariate extension of the concept, with the second moment of the random vector u defined as the matrix E(uu^T), the expected outer product of the activity profiles, where the expectation is across the measurement channels. The second-moment matrix of the activity profiles is given by (2)

Thus, each cell of this matrix contains the scaled inner product of two activity patterns.

Before calculating G, some investigators subtract the mean activity across measurement channels for each condition from the data. In this case, Eq 2 becomes the variance-covariance matrix of the activity profiles –the second moment around the mean activity profile. Here we do not remove the mean, but use the second moment around zero. From the perspective of a neuron that reads out the activity pattern of an area, any difference between activity patterns across conditions can be used to decode information. Some features (for example, stimulus intensity) may be encoded in the mean activity over all measurement channels. Other properties (for example, stimulus identity) may be encoded in relative activity differences, with some measurement channels responding more to one condition, and others to a different condition. The second moment around zero captures both of these potentially meaningful differences.

Any feature of the conditions that we might want to decode can be defined by a K × 1 vector f with one entry per condition, which describes how the feature varies across conditions. To obtain a linear read-out estimate for the feature f_i for a given condition i, we weight each channel’s observed activity using the P × 1 read-out vector v: (3)

We would like the estimate to have very different values for two trials that differ on the feature value, while showing small differences for trials that have the same feature value. f^T U is the pattern that encodes the feature. We are looking for the readout vector v that maximizes the ratio S between the sum-of-squares of the readout of the feature and the sum-of-square of the readout of the noise: (4)

The solution to this equation is commonly known as Fisher’s linear discriminant [35], which, under the assumption of homoscedastic Gaussian noise, is the best achievable linear decoder. If the noise is isotropic (or the data is adequately pre-whitened), then E^T ff^T E = Ib, where b is a constant. The denominator then depends only on the norm of the read-out vector v, not on its direction, and can be ignored when v is constrained to have a norm of 1. The best readout vector v is then given by the first eigenvector of the matrix U^T ff^T U, and the quality of the best readout is determined by the corresponding eigenvalue.

Non-zero eigenvalues (eig) of a square matrix are invariant to cyclic permutations of the product order: (5)

Therefore, the quality of the best linear decoder for any feature (as defined by f) is fully characterized by f^T Gf.

By the same logic, the second moment determines decodeability in general. Assume you want to use a non-linear decoder to estimate feature f. For this to be possible, one has to ensure that one can decode the difference between any two conditions that have different feature values. Once that distinction is made, one can use a arbitrary non-linear function (i.e. table lookup) to read out the feature for each condition. Because the second moment defines the decodeability of any pair of two experimental conditions, it also defines whether a feature can be read out non-linearly.

Importantly, the second moment is only a sufficient statistic for decodeability under the assumption that the readout neurons can integrate information from the entire population that constitutes the code, i.e., it as capable of any arbitrary linear transform of the input data. If the readout neuron is only partially connected, it becomes important to what extent particular information is concentrated in restricted sets of neurons. This information is captured in the higher moments of the activity profile distribution, a point to which we will return in the Discussion. For a fully connected readout neuron that can weight activities in any arbitrary way, the second moment is a sufficient statistic of the decodable information.

Representational analysis in the context of fMRI

The methods in this paper were first developed in the context of fMRI data analysis, and our examples will come from this domain. A simple way to apply the analyses to fMRI data is to use as activity estimates () the regression coefficients, or “beta”-weights, from a first-level time series analysis [36, 37]. The time-series model accounts for the hemodynamic lag and the temporal autocorrelation of the noise. The activity estimates usually express the difference in activity during a condition relative to rest. Activity estimates commonly co-vary together across fMRI imaging runs, because all activity estimates within a partition are measured relative to the same resting baseline. This positive correlation can be reduced by subtracting, within each partition, the mean activity pattern (across conditions) from each activity pattern. This makes the mean of each measurement channel (across condition) zero and thus centers the ensemble of points in activity-pattern space that is centered on the origin.

Rather than using the concatenated activity estimates from different partitions, encoding analysis and PCM can also be applied directly to time series data. As a universal notation that encompasses both situations, we can use a standard linear mixed model [38]: (6) where Y is an N × P matrix of all activity measurements, Z the N × K design matrix, which relates the activity measurements to the K experimental conditions, and X is a second design matrix for nuisance variables. U is the K × P matrix of activity patterns (the random effects), B are the regression coefficients for these nuisance variables (the fixed effects), and E is the matrix of measurement errors. If the data Y are the concatenated activity estimates, the nuisance variables typically only model the mean pattern for each run. If Y consists of time-series data, the nuisance variables typically capture additional effects such as time-series drifts and residual head-motion-related artifacts.

Representational analysis in the context of neurophysiological recordings

All three methods can also be applied to recordings of single cell activity or neurophysiological potentials [9, 25]. The activity estimates can then be firing rates estimated over a temporal window for each trial, or the power in different frequency bands over time. Because the trial-by-trial variability of firing rates will usually increase with the mean firing rate, it is advisable to use the square root of firing rates to make the data conform better to the assumption that the variance of the noise is independent of the signal [39].

Here we focus on models that treat the activity patterns U as static snapshots. To exploit the temporal detail provided by electrophysiological recordings, the analyses can be either performed using a sliding window over the time course of the trial [40–42], or by “stacking” the time series and conditions, resulting in a activity matrix with TK rows [43].

Encoding analysis

An encoding model characterizes the structure of the representation in terms of a set of features [14–17]. We will show in the following that encoding models are representational models as definied by the second moment of the activity profiles. For this to be the case, however, the use of regularized regression is a crticial factor. We will therefore first present the encoding apporach in general, and then show why regularisation is important to test for distributions with a definied second moment.

In general. the value of each feature for each experimental condition is coded in the model matrix M (K conditions by Q features). The feature weight matrix W (Q features by P channels) then determines how the different model features contribute to the activity profiles of different measurement channels to produce the predicted activity patterns U: (7)

Geometrically, we can think of the features as the basis vectors of the subspace, in which the activity profiles reside (Fig 2A).

Encoding analysis without regularization.

To adjudicate among encoding models of different numbers of features—and hence different numbers of free first-level parameters—most researchers use independent test sets [15–17]. A training data set is used to estimate the feature weights for each channel, and the resulting prediction is then evaluated on a held-out test data set. This can be implemented in a statistically efficient manner by using cross-validation, which is usually performed by holding out a single partition (e.g. fMRI imaging run) as a test set, and using the remaining M-1 partitions as the training set. Each partition is held out as the test set once and prediction performance is averaged across the M folds of cross-validation. Encoding models can also make predictions about conditions that are not in the training set (Discussion). However, we focus our simulations on cases, in which training and test sets include the same experimental conditions.

The weights can be chosen to minimize the sum of squared errors on the training data, i.e. using linear regression: (8) where we define to be the average activity estimates from all partitions except m. The prediction for the left-out test data of run m is (9)

The accuracy of the prediction can be assessed by relating the residual sums-of-squares (SSR) of the prediction to the total sums-of-squares (SST) of the observed activities, summed over all partitions, conditions, and voxels (10)

Alternatively, we can evaluate the prediction by correlating the predicted and observed activity patterns across all conditions and channels. Assuming that the mean of each channel across all conditions is zero (given mean pattern subtraction), the correlation is given by (11)

The correlation introduces an arbitrary scaling factor between prediction and observations and, in contrast to Eq 10, allows the model to over- or under-predict the data by a scalar factor without penalty. Encoding analysis can also be applied directly to the time-series data without an intervening model (Eq 6).

To understand how encoding analysis adjudicates between models, consider the graphical representation of the estimation process (Fig 3). In this example, the training data the activity profile of a single measurement channel, which can be visualized as a point in activity-profile space (black cross). Regression analysis can be understood as the orthogonal projection of the measured activity profile onto the linear subspace spanned by the features of the model. The two models depicted in Fig 3A and 3B have different features (blue arrows) that define different subspaces (planes with blue outlines). Therefore, the training data is projected onto two different planes and the prediction for the test data differs between the two models. The model with a subspace that better describes the cloud of activation profiles will make better predictions overall across the measurement channels, show lower cross-validation error, and will hence be more likely selected as the winning model.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Adjudicating between encoding models with and without regularization.

The axes of the three-dimensional space are formed by the response to three experimental conditions. The activity profile of each unit defines a point in this space. Models are defined by their features (blue arrows) and (when using regularization) a prior distribution of the weights for these features. The features and the prior, together, define a distribution of activity profiles (ellipsoids indicate an iso-probability-density contours of the Gaussian distributions). To predict the activity profile of a single measurement channel, the model is fitted to the training data set (cross). Simple regression finds the shortest projection (black dot) onto the subspace defined by the features, whereas regression with regularization (red dot) biases the prediction towards the predicted distribution. Two models (A, B) with features that span different model subspaces are distinguishable using regression without regularization. (C) This model spans the same subspace as model A. Unregularized regression results in the same projection as for model A, whereas regression with regularization leads to a different projection. (D) A saturated model with as many features as conditions. Unregularized regression can perfectly fit any data point (cross and black dot coincide). With regularization, the prediction is biased towards the predicted distribution (iso-probability-density ellipsoid).

https://doi.org/10.1371/journal.pcbi.1005508.g003

Importantly, encoding analysis without regularization compares the subspaces of the competing models, but not their probability distributions. For example, the model depicted in Fig 3C predicts a different distribution than the one in Fig 3A. The features of these two models, however, span the same subspace. Therefore, without regularization, the predictions of these two models are identical (black dots) and the models indistinguishable.

Pattern component modeling

An alternative to cross-validation is to evaluate the likelihood of the measured activity profiles under the representational model. This approach is taken in pattern-component modeling [22]. We start with a generative model of the activity profiles (Eq 6). We consider the activity profiles (columns of U) to come from a multivariate Gaussian distribution with zero mean and second-moment matrix G. To account for other nuisance effects (mean activity for each partition, low-frequency drift, etc), the model also has some fixed-effects regressors (B). We are not interested in fitting U per se, but simply want to evaluate the likelihood of the data under different models, marginalized over all possible values of U. The marginal distribution for each channel (columns of matrix Y) takes the form of a multivariate normal: (16)

The predicted covariance matrix of the activity measurements for each person is the function of the model (as encoded in the second-moment matrix G) and two second-level parameters (θ): one that determines the strength of the signal (s) and one that determines the variance of the noise (). In determining the likelihood, we remove the fixed effects using the residual forming matrix (17)

We need to then account for the removal of these fixed effects by evaluating the restricted likelihood l(Y | G,θ) [47]: (18)

To evaluate the fit of a model, the scaling and noise parameters need to be determined. For fMRI data, these two parameters can vary widely between different brain regions and individuals, and are not meaningful in themselves. We therefore replace θ with point estimates that maximize Eq 18—i.e., the approach uses Empirical Bayes, or Type-II maximum likelihood for model comparison [45]. Because every model has the same two free second-level parameters, even models that are based on different numbers of features can be compared directly. An efficient implementation of this algorithm can be found in the open-source Matlab package for PCM [48].

Representational similarity analysis

Simulated example data sets.

We use simulated data sets here to evaluate and compare the three analysis techniques in a situation where the ground-truth is known. The three simulated example data sets are inspired by real fMRI studies. The first two examples are motivated by a paper investigating the representational structure of finger movements in primary motor and sensory cortex [28]. The structure of the empirically measured distances between movements of the five fingers was highly reliable across different individuals. The main question was whether this invariant structure is best explained by the correlations of finger movements in every-day life—i.e. the natural statistics of movement [50], or by the patterns of muscle activity required for these movements. Rather than hypothesizing that certain features form the basis set generating the activity profiles distribution, we could directly predict the second-moment matrices, and hence the RDMs, from the correlations between naturally occurring movements, or the correlations of muscle activity patterns. The predicted RDM for individuated movements of the five fingers (Exp. 1) is shown in Fig 4A and 4B. The second example comes from experiment 3 in the same paper, this time looking at 31 different finger movements, which span the whole space of possible “piano-chord” combinations (Fig 4C and 4D).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Representational dissimilarity matrices (RDMs) for the models used in simulation.

Each entry of an RDM shows the dissimilarity between the patterns associated with two experimental conditions. RDMs are symmetric about a diagonal of zeros. Note that while zero is meaningfully defined (no difference between conditions), the scaling of the distances is arbitrary. For Experiment 1, the distance between the activity patterns for the five fingers are predicted from the structure of (A) muscle activity and (B) the natural statists of movement. In Experiment 2 (C, D) the same models predict the representational dissimilarities between finger movements for 31 piano-like chords. For Experiment 3 (E, F), model predictions come from the activity of the seven layers of a deep convolutional neural network in response to 96 visual stimuli. The 1^st convolutional layer and the 1^st fully connected layer are shown as examples.

https://doi.org/10.1371/journal.pcbi.1005508.g004

The third example uses an experiment investigating the response of the human inferior temporal cortex to 96 images, including animate and inanimate objects [9]. The model predictions are derived from a convolutional deep neural network model—with each of the 7 layers providing a separate representational model. The bitmap images were presented to the deep neural network and the internal activity patterns used as representational models.

All data sets where simulated with 8 runs, 160 voxels, and independent noise on the observations. The noise variance was set to σ² = 1. We first normalized the model predictions, such that the norm of the predicted squared Euclidean distances was 1. We then derived the second moment matrix (G) for each model using Eq 22 and created true activity patterns that were normally distributed with second moment U U^T/P = Gs. The observation for each run were then generated by adding normally distributed random noise with unit variance to the data (Eq 1). The signal-strength parameter s was varied systematically starting from 0 (pure noise data).

We generated 3,000 data sets for each experiment, parameter setting, and model. Each data set was generated by one model (ground truth) and was analyzed so as to infer the data-generating model, using each of the inference methods. To evaluate how well the methods adjudicated between the models, we compared the fit of the true model (i.e. the model that generated that particular data set) with each alternative model by counting the number of instances, in which the method decided in favor of the correct model. Thus, even though there were 7 alternative models in Experiment 3, chance performance for the pairwise comparisons was always 50%. The percentage of correct decisions over all possible model pairs and simulation was used as a measure of model-selection accuracy.

Encoding analysis without regularization

When evaluating encoding models without using regularization, one compares the subspaces spanned by the respective model features. To make different models distinguishable, one typically needs to reduce the dimensionality of the model matrix M, for example by using only the eigenvectors with the n highest eigenvalues of the predicted second-moment matrix. The decision to use a given number regressors is somewhat arbitrary: For example, Leo et al. [21] used 5 “synergies” (i.e. principal components of the kinematic data of 20 movements), as these explained 90% of the variance of the behavioral data.

Here we explore systematically how the number of principal components influences model selection. For each experiment, we simulated data sets with a fixed signal-to-noise ratio (Exp. 1 and Exp. 3: s = 0.3, Exp 2: s = 0.1; = 1), and compared model selection accuracies using a number of principal components ranging between one and the maximum number. We used both cross-validated performance measures, (Eq 10) and r (the correlation between predicted and observed values; Eq 11) to perform model selection.

Fig 5A–5C shows the percentage of correct model selections for Experiments 1-3. Results for encoding analysis without regularization are shown in blue. The dimensionality that differentiated best between competing models was 2, 3, and 5 features, respectively. As more features were included, the number of correct model selections declined. When the number of features was the same as the number of conditions minus 1 (due to the mean subtraction), i.e. the models became saturated, model selection accuracy fell to chance. This is expected, as two saturated models span exactly the same subspace and hence make identical predictions (Fig 3D).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Dependence of encoding model analysis on regularization and the number of included model features.

(A-C) Percent correct model selections using either (solid line) or correlation (dashed line) for encoding models without a prior (blue lines) and with a prior (red line). (D-F) Correlation between predicted and observed patterns. (G-I) Predictive for the encoding models with prior. All values for models without prior are negative, and therefore not visible.

https://doi.org/10.1371/journal.pcbi.1005508.g005

Using correlations as selection criterion led to more accurate decisions than using . Correlations (Fig 5D–5F, blue lines) were generally positive and peaked at a number of features that was slightly higher than the optimal dimensionality for model selection. values for encoding without a prior were all negative (and are therefore not visible), because the approach does not account for the noise in the data and hence leads to predictions that are too extreme—i.e. the approach over-predicts the scale of the data. Correlations are insensitive to this problem as they allow for arbitrary scaling between predicted and observed values.

Encoding approaches with regularization

From a Bayesian perspective, employing regularization (Eq 13) is equivalent to adding a prior to the feature weights. Note that this changes the representational hypotheses tested. For example, the models for Experiment 3, based on the neural network representations, now predicted not only that some weighted combination of the neural network features can account for the data, but more specifically that the distribution of activity profiles should match the distribution of activity profiles of the original neural network simulation. In the model matrix, we scaled each principal component of G with the square root of the eigenvalue (Eq 15), such that we could employ ridge regression to obtain the best linear unbiased predictor for the held-out data patterns.

For encoding models with a prior, model selection performance increased with increasing number of features (red lines, Fig 5A–5C). Thus, dimensionality reduction of the model is not necessary here. Furthermore, model selection was always more powerful with than without a prior when correlation was used for model selection. This reflects the fact that the prior provides additional information about the models to be compared. It enables us to compare well-defined distributions of activity profiles instead of just subspaces.

For Experiments 2 and 3, the criterion performed substantially worse than the correlation between predicted and observed activity patterns. The difference between the two criteria arises from the fact that correlations allow for an arbitrary scaling between predicted and observed activity patterns, whereas penalizes deviation in scale. The scaling of the prediction in turn strongly depends on the choice of the scalar regularization coefficient. This fact is illustrated in Fig 6, where we simulated data from Exp. 2 with a fixed noise and signal strength, and varied the regularization coefficient systematically. While is highly sensitive to the choice of the regularization coefficient, the correlation criterion is not. Because the regularization coefficient is determined separately for each cross-validation fold and model, deviations from the optimal ridge will decrease model selection accuracy for the cross-validated criterion, but not for the correlation criterion.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Sensitivity of the (solid line) and correlation (dashed line) to the choice of the regularization coefficient.

Simulations come from Experiment 2 with a true signal strength of s = 0.2 and a noise variance of 1. For this combination the optimal regularization coefficient is (dashed vertical line). The correlation criterion is generally robust against non-optimal choice of regularization coefficient.

https://doi.org/10.1371/journal.pcbi.1005508.g006

In sum, using regularization improves model selection performance, even if the encoding model has fewer features than conditions or measurements. Rather than just comparing subspaces, the implicit prior on the weights means that a more specific hypothesis is being tested. From this perspective, it is unsurprising that we can adjudicate between these hypotheses with greater accuracy. Furthermore, the use of correlation instead of the predictive makes model selection more robust against variations in the regularization coefficient.

Representational similarity analysis

When evaluating models with RSA, there is no need to restrict the model to a specific number of features—the second-moment matrix from all features can determine the predicted distances. As an empirical dissimilarity measure, we used the crossnobis estimator [32] and compared the predicted to the measured RDM. To select the winning model, we used rank-based correlation of dissimilarities [27], Pearson correlation, correlation with a fixed intercept (Eq 24), and the likelihood of the observed distances under the normal approximation (Eq 26) using the full variance-covariance matrix of the estimated dissimilarities.

For Experiment 1 (Fig 7), rank-based correlation performed substantially worse than the other criteria. The lower performance of rank correlation may have been exacerbated here by the fact that the two models predict relatively similar dissimilarity ranks. However, we expect lower performance for rank correlation in general, because this approach does not use all the information in the measured RDMs. It forgoes the assumption of a linear relationship between predicted and measured dissimilarities and therefore does not exploit the information in the continuous magnitudes of the dissimilarities. Likelihood-based RSA yielded the best decisions; slightly better than Pearson correlation and fixed-intercept correlation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. RSA model selection accuracies for different criteria of RDM fit.

Data sets for all three experiments were generated with varying signal strength (horizontal axis). The percentage of correct decisions using different criteria is shown (dotted line). Models were selected based on the Spearman rank correlation (purple), Pearson correlation (green), fixed intercept correlation (blue) or likelihood under the multinormal approximation (red). For comparison, the model selection accuracy for PCM is shown in the dotted line.

https://doi.org/10.1371/journal.pcbi.1005508.g007

The advantage of the likelihood-based approach was clearer for Exp. 2 and 3. Here, it led to about 10 percentage points greater accuracy of the decisions than the next-best RSA approach. This advantage is likely due to the fact that Pearson correlations and especially fixed-intercept correlations (Eq 24) are sensitive to the observed value for the largest predicted dissimilarities, as these data points have a large leverage on the estimated regression line. Indeed, some of the models for Exp. 2 and 3 contain a few especially large dissimilarities, which will influence the model fit strongly. The likelihood-based approach incorporates the knowledge that large dissimilarities are measured with substantially larger variability [29], and hence discounts their influence. Notably, rank-based correlation performed relatively well on these models as compared to Pearson correlation, likely because rank correlation is robust to outliers and less dominated by the large predicted distances.

In sum, these simulations show that it is advantageous to take the covariance structure of the measured dissimilarities into account whenever the additional assumptions this requires are justified.

Pattern component modeling

In the same simulations, we also applied the direct likelihood-ratio test, as implemented by PCM. As all the assumptions of the generative model behind PCM are met in the simulation, we would expect, by the Neyman-Pearson lemma [24], that this method should provide us with highest achievable model selection accuracy. Model selection performance (dotted line in Fig 7) was indeed systematically higher than for the best RSA-based method. For direct comparison of the so far best methods—PCM, likelihood-based RSA, and encoding analysis with regularization (using correlations as a model selection criterion)—we simulated the three Experiments at a single signal strength (Fig 8).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Model selection accuracy and execution time for likelihood-based RSA, PCM, and encoding analysis with regularization.

(A-C) Model-selection accuracy was inferentially compared between the three techniques on the basis of N = 3,000 simulations, using a likelihood-ratio test of counts of correct model decisions [51]. The signal-strength parameter for the simulation was set to s = 0.3 for Exp. 1, s = 0.15 for Exp. 2, and s = 0.5 for Exp. 3. All resulting significant differences (two-tailed, p<0.01, uncorrected) are indicated by a horizontal line above the bars. (D-F) Execution times for the evaluation of a single data set under a single model. For encoding, the time is split into the time required to estimate regression coefficients (dark blue) and the time to determine the regularization constant (light blue).

https://doi.org/10.1371/journal.pcbi.1005508.g008

In this simulation, PCM resulted in 1.48, 3.01 and 2.86 percentage points (for Exp. 1-3, respectively) better model selection accuracy than likelihood-based RSA, and 1.98, 1.17 and 0.85 percentage points higher model selection accuracies than an encoding analysis using correlations. PCM never performed worse than another method and performed significantly better than the other two approaches in 4 of 6 total comparisons across the three experiments (Fig 8). There were no significant performance differences between RSA and encoding analysis. Overall, however, all three methods were very close in performance.

Computational cost

A practical concern is the speed at which the model comparison can be performed. This is usually not important when evaluating the model fit on a small number of participants or ROIs. However, if a larger number of models is evaluated continuously over the cortical surface using a searchlight approach [52, 53], or in data sets with large numbers of participants, computational cost becomes a practical issue. While we cannot treat this issue exhaustively, we provide here a brief overview over the computation time required for the three methods for our specific examples and implementation. In general, the computation time will of course depend on the number of conditions, the number of channels, the exact variant of each technique. Our goal here is simply to give the reader a starting point for making a choice for a particular application, trading off computational and statistical efficiency.

Both RSA and PCM operate on the inner product matrix of the activity estimates, thus the computational costs for these approaches is virtually independent of the number of voxels. PCM works on the MK × MK inner product matrix of the activity estimates, whereas RSA operates on a K × K matrix of distances between conditions. For a small number of conditions, this explains the favorable computational cost of RSA. However, when using likelihood-based RSA, the covariance matrix of the distances needs to be calculated and inverted. The size of this matrix is (K(K-1)/2)² and it therefore grows with the 4^th power of the number of conditions K. For Exp. 3 (Fig 8F) with K = 96, this is computationally costly, whereas PCM still only needs to invert matrices of size (MK)². Using RDM-correlation-based model selection (whether rank, Pearson or fixed-intercept), RSA is much more computationally efficient (not shown).

For encoding models, conducting the actual ridge regression for each cross-validation fold (dark blue area) is extremely fast and efficient. The main cost of the technique lies in the determination of the optimal ridge coefficient (light blue area). In our simulations, we use restricted maximum likelihood estimation (Eq 18) to do so—therefore this cost is always M times higher than for PCM alone. Depending on the implementation, generalized cross-validation [46] may offer a considerable speed-up. If very high speeds are required, one could use a constant ridge coefficient and accept the possible loss in model selection accuracy. In sum, while PCM is computationally feasible across the three experiments, encoding models were less efficient in the present implementation and likelihood-based RSA was less efficient than PCM for the condition-rich scenario of Experiment 3. Alternative variants of encoding models (with fixed ridge coefficient) and RSA (with correlation-based model selection) are less statistically efficient, but beat PCM in terms of computational efficiency.

Encoding models without prior define subspaces, not distributions of activity profiles

There is a fundamental difference between encoding models with and without weight priors. Without a prior on the feature weights, encoding models test how well the subspace spanned by the model features captures the observed activity profiles. For models to be discriminable, the dimensionality (i.e. the number of features) of each model must be substantially lower than the number of experimental conditions. As the number of model dimensions increases, the subspaces of competing models increasingly overlap. Once the number of features matches the number of experimental conditions, their subspaces comprise the entire space of activity profiles, each perfectly fits the training data, and their predictions for unseen data become identical.

A subspace specifies what activity profiles are possible and what activity profiles are impossible (though they might still arise as estimates because of the noise). A subspace might be conceptualized as an infinite flat distribution over the subspace dimensions, with 0 probability outside the subspace. However, a uniform distribution on an infinite interval has an infinite second moment and hence does not specify the neural representation uniquely.

L2-norm regularization (i.e. ridge regression) is equivalent to imposing a Gaussian prior on the regression weights. With such a prior, the representational model specifies a probability distribution with a finite second moment. When changing the form of regularization, one also changes the implicit prior, and hence the representational model that is being tested. Thus, regularization is not simply a trick for stabilizing the fit. Instead, the weight prior forms an integral part of the model, which determines the strength with which each feature is encoded according to the model. Choosing a specific form of regularisation therefore constitutes a decision about the neuroscientific hypothesis to be tested rather than a methodological consideration.

Encoding models tests hypotheses about activity profile distributions, not features sets

Encoding models do not support inferences about the particular feature set generating a representation, because infinitely many feature sets can span the same space. Even when using a prior, the feature set that characterizes a given representational model is not unique. Features should not in general be constrained to be orthogonal in the space of experimental conditions, because the structure of the model is not usually meant to depend on the experimental conditions chosen. Whether the features chosen are orthogonal or not, there is an infinite number of basis sets of features that express the same representational model (inducing the same second moment of activity profiles, Eq 3). For example, two equally long correlated feature vectors can equally well describe a distribution with elliptical isoprobability-density contours (Fig 3A) as two orthogonal features, with one vector longer than the other. Thus, when one representational model is shown to be superior to others, it does not imply anything special about the feature set chosen to express that model. These complications need to be kept in mind in the interpretation of the results of encoding model analyses. It is very tempting to attribute meaning to the particular features, especially when they are mapped onto the cortical surface [17, 21]. When interpreting these maps, one needs to remember that a feature set only describes a distribution of activity profiles, and that very different maps can emerge when the same distribution is described by a rotated set. In PCM and RSA, the equivalence of different feature sets is made explicit, as they lead to the same second-moment and representational dissimilarity matrices.

Likelihood-based RSA is more sensitive than correlation-based RSA

When using RSA to test representational models, the crossnobis estimator provides a highly reliable measure of dissimilarity with the added advantage of having an interpretable zero-point [32]. Rank-based, Pearson, and fixed-intercept correlation provide fast and straightforward ways of measuring the correspondence between predicted and observed distances, so as to select the representational model most consistent with the data. However, using simple correlations ignores the dependence of the distance estimates, as well as their unequal variances. In other words, the sampling distribution of the estimated RDM in the space spanned by the dissimilarities (one dimension per pair of conditions) is not isotropic. This problem is addressed in likelihood-based RSA, which uses a multivariate-normal approximation to the sampling distribution of the crossnobis RDM estimate [29]. The approximation provides an analytical expression for the statistical dependency of distance estimates, as well as their signal-dependent variances. In the simulations, likelihood-based RSA was shown to be more powerful than correlation-based RSA. Its model-selection accuracy was only slightly below the theoretical upper bound, as established by PCM. Likelihood-based RSA might therefore become the approach of choice when comparing representational models using crossnobis estimates.

There are situations, however, in which the models are not specific enough to support ratio-scale predictions of representational dissimilarities. Moreover, for measurement modalities like fMRI, it might be undesirable to assume a linear relationship between predicted and measured representational dissimilarities. Rank-correlation-based RSA [25, 27] provides a robust method that is not dependent on the assumption of a linear reflection of the underlying neural dissimilarities in the data RDM. It is also more computationally efficient in the context of condition-rich designs. Likelihood-based RSA becomes computationally expensive as the number of conditions increases. A practical compromise might be to only use the diagonal of the variance-covariance matrix, which would dramatically reduce computational complexity at the expense of neglecting dependencies among dissimilarity estimates.

Which method is best?

For all simulations, model selection using PCM [22] was better than competing methods. This is not surprising, as the data were simulated exactly according to the generative model underlying this approach (Gaussian distribution of noise and signal, independence across voxels). In this case, PCM implements the likelihood-ratio test, which by the Neyman-Pearson lemma [24] is the most powerful test. Beyond confirming what we know from theory, the simulations were important because they revealed how close the other two techniques come to the theoretical upper bound established by PCM. Results showed that encoding models with a prior and likelihood-based RSA perform near-optimally. In practice, we therefore expect these three techniques to provide similar answers. When its assumptions hold, PCM has clear advantages for model comparison, providing optimal power at reasonable computational cost. However, the other two techniques have other advantages that make them attractive for specific applications.

RSA using RDM correlation for model selection gives up statistical efficiency for computational efficiency, and beats PCM at the latter. When rank correlation is used to compare RDMs, the inference does not rely on a linear relationship between the true dissimilarities and the estimated dissimilarities, an assumption that might be violated in many contexts. RSA also provides readily interpretable intermediate statistics (cross-validated distances), which are closely related to linear decoders for all pairs of stimuli. These statistics can be used to test whether two conditions have different activity patterns [27, 29], or whether the dissimilarity is larger for one pair than for another pair of conditions. Multidimensional scaling of the stimuli on the basis of their representational dissimilarities also provides an intuitive visualization of the representational structure [25], which can be very helpful in the generation of novel representational hypotheses.

In contrast, PCM sometimes demands complicated approaches to answer simple questions: For example, to test the hypothesis that a difference between two conditions is encoded, one would need to fit one model that allows for separate patterns and one model that does not—and then compare the marginal likelihood of these models. Furthermore, PCM requires the noise to be explicitly modeled, whereas RSA removes the bias arising from noise through cross-validated distances.

Encoding analysis explicitly estimates the first-level parameters that describe the response for each individual voxel. This enables the mapping of the estimated features onto the cortical surface to study their spatial distribution [17, 21].

In sum, the three methods are deeply related in that they test hypotheses about the second moment of the activity profiles. However, each method constitutes a unique perspective on the data and supports different kinds of exploratory analyses. We view the methods as complementary tools that are part of a single coherent toolkit for analyzing representations.

Single-voxel vs. multi-voxel inference

An important issue, which we have not touched upon so far, is whether to perform model comparison on single or multiple voxels. While RSA and PCM are usually applied to groups of voxels (such as for ROIs or searchlights), encoding models are often compared on the single-voxel level. This tendency, however, is not strictly inherent in methodological constraints: The searchlight approach for RSA and PCM can be reduced to single voxels, and encoding models can be combined with multi-voxel searchlights. Analyses with coarser granularity give up some spatial precision of the map in exchange for greater statistical power. Searchlight mapping boosts power (1) by locally combining the evidence, (2) by enabling the use of a multivariate noise normalization, and (3) by reducing the effective number of multiple comparisons [54]. There is no reason to assume that a single-voxel searchlight is always the optimal choice when balancing spatial precision and power. Based on our previous results [32], we expect that ignoring voxel dependencies will entail a loss of sensitivity when making inferences on representational models for regions of interest comprising multiple responses.

Testing models without overfitting to the noise and to the sample of experimental conditions

Whenever a model is fitted using experimental data, its parameters will necessarily be overfitted to the data to some extent. Assessing the performance of a fitted model therefore requires independent test data. An important question is whether the test data should consist in independent measurements for the same experimental conditions or in measurements for a fresh sample of experimental conditions (e.g. a different sample of visual images). The simple answer is that it depends on the inference we would like to make. If our hypothesis is restricted to the present set of conditions (e.g. five finger movements), we need only account for overfitting to the noise in the data and require different measurements for the same conditions. If our hypothesis is about a population of conditions (e.g. all natural images), we need to account for overfitting to the condition sample and require measurements for an independent random sample of conditions from the population of conditions covered by our hypothesis.

However, overfitting only needs to be accounted for when the model being tested had parameters fitted in the first place. Encoding models always require independent test sets to account for the over-fitting of the first-level parameters of the representational model (feature weights). RSA and PCM, by contrast, rely explicitly on summary statistics of the responses. Therefore, only second-level parameters related to the strength of the signal and noise need to be fitted (see Table 1). Because the representational models considered here had the same number of such second- level parameters, they could be compared directly.

What about decoding approaches?

Decoding is widely used in multivariate analysis of brain imaging data [11–13]. Can it serve us also as a tool for comparing representational models? While one can use standard decoding approaches to determine whether specific features are represented in an area or not, it does not lend itself to the comparison of full representational models (as defined here). Representational models determine (via the second moment matrix) the decodability of any linear feature, not just a restricted set of features. This is most obvious in RSA, where the RDM assembles all pairwise condition discriminabilities. It is of course possible to use decoding in the context of the methods considered here. For example, some studies have used encoding models to decode stimuli [15, 16, 21]. Decoding accuracy on held-out data can then serve, instead of correlation or , as a performance evaluation of an encoding model. While this approach is motivated by the intuitive demonstration of mind reading, it does not provide a particularly natural or powerful approach to adjudicating between representational models. Alternatively, we could use classification accuracy as a measure of dissimilarity between two conditions in the context of RSA [55]. However, classification essentially converts a continuous measure of dissimilarity into a binary label of correct / incorrect. It is therefore expected to be less informative than the underlying continuous measure, and we have shown previously that this entails a loss of sensitivity in practice [32]. In sum, decoding is not particularly useful for the evaluation of representational models [14, 23] and should therefore be limited to situations, in which the quality of the decoding itself is the measure of interest.

Flexible representational models

All models considered here were “fixed”, i.e., they did not include free parameters that would change the predicted second-moment matrix. In many applications, however, the relative importance of different features (for example encoding strength for orientation and color) are unknown. In this case, the predicted second moment can be expressed as the weighted sum of different pattern components, i.e. G = ∑_i ω_i C_i [22, 56–58], with the weights being free second-level parameters. In other situations, G is a nonlinear function of free model parameters: For example, G depends non-linearly on the spatial tuning width in population receptive field modeling [59]. Both RSA and PCM already provide a mechanism to estimate such parameters, as both approaches need to estimate the signal strength parameters s by maximizing the respective likelihood function (Eqs 17 and 28)—and the analytical derivatives of the likelihood (Eqs 17 and 28) with respect to the parameters are easily obtained. In the context of encoding approaches using ridge regression, free model parameters that change the model structure would result in independent scaling of different features, rotations, or extensions of the model matrix M. At the time of writing there are no published examples of such parameter optimization in the context of cross-validated encoding models that we know of.

The inclusion of free parameters into the model also enables the specification of measurement models. Representational models ideally test hypotheses about the distribution of activation profiles of the core computational elements—i.e. neurons. When using indirect measures of brain activity such as fMRI or MEG, the distribution of activity profiles across measurement channels is also influenced by the measurement process, which samples and mixes neuronal activity signals in the measurement channels [30, 60–63]. As the underlying brain computational models become more specific and detailed, the corresponding measurement models will also have to be improved.

Higher-order moments of the activity profiles

We focused on approaches that characterize the distribution of activity profiles by its second moment. If the true distribution of the activity profiles is a multivariate Gaussian, then the second moment fully defines the distribution of activity profiles. However, a representational hypothesis may not only predict that the response for condition A is uncorrelated to the response for condition B, but, for example, that channels either respond to A or B, but not to both A and B. Such tuning is for example prevalent in primary visual cortex, where neurons (and voxels) respond to a stimulus in a one specific part of the visual field, but less often two or more disparate locations [59]. This would correspond to a non-Gaussian prior on the feature weights. In a recent publication, Norman-Haignere and colleagues [64] suggested a likelihood-based method, in which the Gaussian prior on the feature weights W is replaced with a Gamma distribution, essentially providing a non-Gaussian extension of PCM. It will be interesting to determine to what degree such non-Gaussian distributions are present in fMRI or single-cell data, and what computational function these may serve.

It is important to stress that the approaches considered here are still appropriate when the distribution of activity profiles is truly non-Gaussian. Even in the non-Gaussian case, the second moment determines the representational geometry and thus the decodability of all possible features. It therefore remains essential for characterizing the representation. Taking into account higher moments of the activity profile distribution would enable us to distinguish between representations that afford the same decoding of features (assuming that readout neurons have access to the entire code), but achieve this by distinct population codes.

Conclusions

If advances in brain-activity measurements are to yield theoretical insights into brain computation, they need to be complemented by analytical methods to test computational models of information processing [65]. The main purpose of this paper was to provide a clear definition of one important class of models—representational models—and to compare three important approaches of testing these. We have shown that PCM, RSA and encoding analysis are all closely related, testing hypotheses about the distribution of activity profiles. Moreover, all three approaches, in their dominant implementations, are sensitive only to distinctions between representations that are reflected in the second moment of the activity profiles. Thus, these three methods are properly understood as components of a single analytical framework. Each of the three methods has particular advantages and disadvantages and preferred areas of application.

1. PCM provides an analytic expression for the marginal likelihood of the data under the model, and therefore constitutes the most powerful test for adjudicating between representational models if the assumptions hold. Its analytical tractability and relative computational efficiency are further attractive features, especially when considering models with increasing numbers of free parameters.
2. RSA provides highly interpretable intermediate statistics and is therefore ideally suited for the visualization and exploratory analysis. Furthermore, simple models are often more easily tested than with PCM. The normal approximation to the distribution of estimated distances enables inference that is nearly as powerful as the likelihood-ratio test provided by PCM. Finally, dissimilarity-rank-based RSA, though less sensitive, provides a means of inference that does not rely on the assumption of a linear relationship between predicted and measured dissimilarities and is computationally efficient even for condition-rich designs.
3. Encoding approaches enable the voxel-wise mapping of model features onto the cortical surface. They therefore are the natural choice when the spatial distribution of features or the voxel-wise comparison of representational models is the main interest.

We hope that the general framework presented here will enable researchers to combine these approaches to make progress revealing the computational mechanisms of biological brains.