Overcoming bias in representational similarity analysis
Roberto Viviani
1: Institute of Psychology, University of Innsbruck, Austria 2: Psychiatry and Psychotherapy Clinic III, University of Ulm, Germany Roberto Viviani University of Innsbruck Institute of Psychology Innrain 52 6020 Innsbruck, Austria [email protected]
Abstract
Representational similarity analysis (RSA) is a multivariate technique to investigate cortical representations of objects or constructs. While avoiding ill-posed matrix inversions that plague multivariate approaches in the presence of many outcome variables, it suffers from the confound arising from the non-orthogonality of the design matrix. Here, a partial correlation approach will be explored to adjust for this source of bias by partialling out this confound. A formal analysis will show the dependence of this confound on the temporal correlation model of the sequential observations, motivating a data-driven approach that avoids the problem of misspecification of this model. However, where the autocorrelation locally diverges from the volume average, bias may be difficult to control for exactly (local bias), given the difficulties of estimating the precise form of the confound at each voxel. Application to real data shows the effectiveness of the partial correlation approach, suggesting the impact of local bias to be minor. However, where the control for bias locally fails, possible spurious associations with the similarity matrix of the stimuli may emerge. This limitation may be intrinsic to RSA applied to non-orthogonal designs. Overcoming bias in representational similarity analysis
Introduction
Representational similarity analysis (RSA, Kriegeskorte et al. 2008) is a multivariate technique to investigate cortical representations of objects or constructs defined by a set of properties (the stimuli). A defining characteristic of this approach is that it eschews inversion of estimates of large covariance matrices, an ill-posed problem that easily gives rise to unstable estimates. Instead, RSA seeks evidence of cortical encoding of properties of stimuli by assessing the concordance between the matrix encoding the pairwise similarity of stimuli and a matrix encoding the similarity of the brain responses to the same stimuli (in an equivalent alternative formulation, the concordance is based on measures of dissimilarity). Combined with the searchlight approach (Kreigeskorte et al. 2006), RSA allows mapping the cortex to seek evidence for the representation of stimuli based on their properties as defined by the experimenter. A downside of RSA is that, unless the coefficients of the model of the brain responses were obtained with an orthogonal design, the covariance of these coefficients is confounded by the covariance imposed by the design. Cai et al. (2019) have drawn attention to this problem, and proposed a maximum-likelihood/Bayesian approach to overcome it. As in another related approach (Friston et al. 2019), this requires knowing the form of the confound to model it parametrically. In contrast with this literature, it will not be assumed here that the form of the confound is known, instead exploring the effectiveness of estimating it directly from the data of the whole volume. This approach is reminiscent of PET techniques of adjusting model coefficients for global signal levels, except that the global term here is a covariance, and adjustment takes place via partial correlation. This paper is structured as follows. The composition of the sum of squares and cross-products matrix will be introduced in the simple setting of multivariate ordinary least squares. Using simulated data, the existence of bias will be demonstrated and the partial correlation approach to overcome this bias. We will then consider generalized least squares models and raise the issue of estimating this confound in a realistic setting. A key issue here will be adequately modeling the autocorrelation of errors in the linear model, i.e. the lack of independence of sequential observations. The proposed approach will be applied to two samples of real EPI data, one collected at ordinary and the other at short TRs (where lack of independence of sequential observations may be pronounced). The Discussion will mention the limitations of the RSA with the partial correlation approach. The Material and Methods section will contain details on the acquisition parameters used in the MRI data. The approach presented here was developed as an add-on function to the SPM package that will be made publicly available. Results
Theory and simulations, ordinary least squares
When collected from a searchlight volume, the model may be written
Y XB E , where Y is the n scans p searchlight voxels matrix of the imaging data, X is a n q design matrix usually containing q BOLD-convolved regressors for the signal elicited by the presentation of stimuli during the experiment, B is a q p matrix of coefficients to be estimated by least squares, and E the n p matrix of errors. This model is applied repeatedly to each searchlight; to simplify notation, indexes referring to the searchlight are omitted. It will also be initially assumed that the rows of E be independent, an assumption that will be relaxed later. Well-known results of multivariate regression are that the least squares estimates of B are identical to those obtained by regressing on the individual columns of Y separately (as in "the massively univariate approach" commonly adopted in functional imaging analyses), and that the expectation of these estimates and the expected variance-covariance are given by ˆ( ) V V E B B , ˆVar( ) ( ) V B Σ X X . (1) In this expression, ( ) V is the operator that stacks the columns of a matrix upon each other into a vector, and is the variance-covariance of the columns of the errors (i.e., over voxels). From equation (1), the expected variance of a column ˆ i β of ˆ B (the coefficient estimates from voxel i) is ˆVar( ) ( ) i i β X X , (2) where i is the diagonal term of giving the variance of the errors in voxel i (Mardia et al. 1979, pp. 160-161, 180-181). After estimating the coefficients ˆ i β , the RSA used here proceeds by forming their sum of squares and cross-products matrix ˆ ˆˆ ˆ i i i BB β β before assessing the correlation of its off-diagonal terms with the analogous terms of the similarity matrix of the stimuli. There are a few different approaches in the literature at this step. The rsa toolbox developed as the original implementation of this technique (Nili et al. 2014) centers the coefficients to compute a dissimilarity matrix based on their pairwise correlations. The sum of squares and cross-products is used here because its involvement in RSA bias is most direct in this case, facilitating its correction. Note also that ˆ ˆ BB is not a sum of squares and cross-products matrix in the usual sense, because the individual ˆ i β are not independent from each other (as specified by in eq. 1), being measured from nearby voxels. Based on the regression on the individual columns of Y, another relevant standard result is the following,
11 1 ˆ ( )( ) ( )( ) i i i ii i β X X X yX X X Xβ εβ X X X ε where i ε is a column of E (Mardia et al. 1979, pp. 160). Taking the squares of all terms and summing over i, one obtains ˆ ˆ ( )ˆˆ ˆ( ) Var( ) i i i i ii i i ii E β β β β X XBB BB β (3) showing that, as an estimate of BB , the off-diagonal terms of ˆ ˆ BB are confounded by those of ˆVar( ) i β , unless the columns of X are orthogonal and ( ) X X I . RSA bias arises from an accidental correlation between the similarity matrix of the stimuli and ( ) X X . Note that the confound due to ˆVar( ) i β is the same over all voxels and searchlights up to a multiplicative factor given by the voxel error variances i . Since also the similarity matrix of the stimuli is the same over all searchlights, it gives rise to a global bias whose intensity is modulated by the voxel error variances but has otherwise the same origin over the whole volume. In contrast, BB will usually vary between searchlights, and its off-diagonal terms will be correlated with the similarity matrix of stimuli only in few sparse areas. This suggests that the average RSA correlation from the searchlights over the whole volume may provide a rough but simple diagnostic for the existence of global bias.
0 45 90 136time 0 45 90 136time pattern A pattern B pattern A pattern B pattern A pattern B−0.100.10.20.3 M ean c o rr e l a t i on i n v o l u m e RSA RSA partial correlation
A B C
Figure 1. A: Schematic illustration of the design used to elicit bias in RSA. Classification of trials in two categories (shown in red and blue) was chosen to produce correlation of coefficient models between trials of different categories (‘pattern A’, left) or between trials of the same category (‘pattern B’, right). B: representational similarity analysis for patterns A and B in volumes containing white noise, boxplots of estimated correlations. On the y axis the average correlations over the whole volume in standard RSA (left) and for partial correlations with the off-diagonal term of (X'X) -1 as covariate. s.d.: standard deviation. The global bias arising from the non-orthogonality of the regressors is illustrated in Figure 1, assessed by the average correlation over the volume in simulated noise data with parameters comparable to those of a realistic MRI dataset and a design with strong non-orthgonality (N = 30; for details for the parameters used in the simulation, see the Methods section). In panel A, BOLD-convolved regressors for trials modeling exposure to stimuli for three seconds are shown, grouped into 4 blocks of 4 trials each. Pattern A, on the left, comprises stimuli of two different categories (in red and blue) shown in alternate series within the blocks. In pattern B, the stimuli within the blocks belong to the same category, with alternate blocks. While this design may not be representative of those found in typical RSA studies, it provides a good challenge for the approaches to correct for the bias that follow. Figure 1B shows boxplots of the mean volume correlation obtained by conducting a searchlight RSA of the similarity matrices of these two patterns after fitting the model to white noise data, where BB is zero. One can see that pattern A resulted in artificially inflated positive correlations, while pattern B produced the opposite effect. Bias was apparent also when using the original rsa toolbox (Nili et al. 2014), where we obtained a mean correlation of 0.19 (std. dev. 0.0043) for pattern A and -0.10 (std. dev. 0.0003) for pattern B. Using a permutation strategy by randomly varying the assignment of categories to trials to compute significance (Kriegeskorte et al. 2008) will be of little help here, since patterns A and B represent cases in which categories are ordered so as to capture the patterns induced by the design matrix best, albeit in different directions. Most permutations will deliver category patterns that are some mixture of patterns A and B, resulting in global bias values intermediate between those obtained with these two patterns. Hence, the correlations obtained with patterns A and B will be located at the far tails of the permutation distribution even if they were computed from noise. To remedy the effects of the covariance of regression coefficients on RSA, the approach explored here consists in partialling out the effects of ( ) X X on the off-diagonal terms of ˆ ˆ BB using a partial correlation, instead of using a simple correlation to test for the presence of the similarity pattern in ˆ ˆ BB . Figure 1C shows that this approach succeeds in taking the average volume RSA correlations to zero in the white noise data. In this simple case, this is almost a given, because the confound is exactly known up to the multiplicative factor of voxel error variance. Generalized least squares
Recall that these results were obtained under the simplifying assumption that the observations were independent. Most functional imaging analysis approaches reject this assumption and fit a generalized least squares model, for example by including a first-order autoregressive term G estimated from residuals pooled over voxels, ˆ ( ) i i β X G X X G y . The exact form of G depends on the model for the temporal autocorrelation of residuals and on implementation details (in the SPM package used in the present work, G V H V , where V is an estimated AR(1) model for the residuals and H a residual-forming projection matrix implementing a high-pass filter). If this estimate is a good estimate of the true temporal dependency structure of the residuals, a standard result is that ˆVar( ) ( ) i i β X G X , (4) which replaces the corresponding term in eq. (3). In real datasets, where generalized least squares is used, one may then use the off-diagonal terms of Bcov = ( ) X G X to partial out the effects of the confound due to the non-orthogonal design. In practice, however, the autocorrelation model may not capture the true structure of the autocorrelation of the data. Furthermore, its estimate is usually computed from residuals pooled from multiple voxels, ignoring possible variations of the autocorrelation between voxels. Hence, a more realistic expression for the expected variance-covariance of the model coefficients may be given by ˆVar( ) ( ) ( ) i i i β X G X X G Γ G X X G X (5) (Diggle et al. 2002, p. 60), where i is the true dependency of the errors in voxel i. Only if G is a good estimate of i , this expression simplifies to ( ) i X G X . This shows that, unless the autocorrelation term is appropriately modeled, residual confounding from the non-orthogonal design may still be present even after partialling out the effect of the off-diagonal terms of
Bcov . Unfortunately, appropriately modeling error dependencies may be an arduous task because the true structure of i is unknown (i.e, whether AR(1), AR(2) or something else) and may be influenced by study-specific factors (such as the properties of the MRI sequence, of the coil, and the frequencies of the design matrix). As before, most terms of this expression are identical over the volume, inducing global bias, while the voxel-by-voxel variation of i introduces the possibility of local bias. For the partial correlation approach, local bias means that one adjustment term alone may not control for the non-orthogonality of the design matrix satisfactorily. Nevertheless, even if G is not a good approximation of i , it might be possible to address global bias by estimating the average variance-covariance * ˆVar( ) β from the whole volume: * v i ii v β β β β β , v ii v β β . (6) Here, the asterisk in the subscript indicates that the variance-covariance is estimated from the whole volume (as opposed to any individual searchlight), and v is the number of voxels in the volume. One may then use the off-diagonal terms of * ˆVar( ) β Svar in a partial correlation to control for global bias. The rationale for this approach is that global bias may be dominant given that most of the terms of eq. (5) are constant across the whole volume, and may be estimated as the average variance-covariance. A similar, even simpler approach consists of using the off-diagonal terms of ˆ ˆ v B B BB to control for global bias, where the asterisk in the subscript of * ˆ B indicates again that this is the matrix of the coefficients for the whole volume. BB is confounded by the mean * BB over the whole volume (eq. 3). However, since the local BB may be sparse and vary over searchlights, they average out, with BB having an effect equivalent to Svar in controlling from global bias (see the Results section below).
Partial correlation in EPI data
Global bias and its correction with partial correlation RSA was explored in real MRI data using the same design as in Figure 1A (Labek et al. 2017, N = 27, TR = 2.27 sec; for details, see the Methods section). The searchlight approach was used to seek evidence for representations of pattern A and B (Figure 2).
Figure 2. Similarity matrices for pattern A and B (left) and confound terms
Bcov and BB (right) averaged over the whole sample. The data were here fitted with an AR(1) model for the temporal correlation of residuals. The first off-diagonal of Bcov and BB is similar to that of pattern A, while containing mostly values of the opposite sign than pattern B. Note that BB , unlike Bcov , is computed without any knowledge of the design matrix.
In Figure 3, panel A, global bias is assessed as before by looking at average volume correlations. On the left ( ptn wo adjustment ), there is evidence of global bias, similarly to the simulations conducted on noise data.In panel B, one can see second-level t maps of the RSA correlation for pattern A. The underlying average RSA correlations are positive in the whole volume, with peaks reaching r = 0.15. This may not seem much but the t map at the second level gave values ranging between 8 and 25 across the volume, which would lead to rejecting the null hypothesis in all voxels. In panel C, left, one can see the RSA correlations computed with the searchlight approach of Bcov itself, as estimated from the AR(1) model of the residuals. One can see that the spatial distribution of the artefactual effect of pattern A, affecting prevalently the white matter region, is similar to that of
Bcov in panel C. This is the behavior described in eq. (3) and (4), and may be understood by noting the similarity of the first off-diagonal of the similarity matrix of pattern A with the same off-diagonal of
Bcov (Figure 2). Inspection of the artefactual correlation for pattern B revealed a similar effect over white matter, but in the negative direction (not shown here for brevity).
C F
Bcov BB|Bcov Bcov BB|Bcovr
B E pattern A pattern At -20 -10 0 -20 -10 0 -20 -10 0 -20 -10 0 -20 -10 0 -20 -10 0 pattern A pattern B−0.100.10.20.3 M ean c o rr e l a t i on i n v o l u m e ptn wo adjustment pattern A pattern B−0.100.10.20.3 M ean c o rr e l a t i on i n v o l u m e ptn | Bcov pattern A pattern B−0.100.10.20.3 M ean c o rr e l a t i on i n v o l u m e ptn | BB pattern A pattern B−0.100.10.20.3 M ean c o rr e l a t i on i n v o l u m e ptn | Bcov+BB FAST D pattern A pattern B−0.100.10.20.3 M ean c o rr e l a t i on i n v o l u m e ptn wo adjustment pattern A pattern B−0.100.10.20.3 M ean c o rr e l a t i on i n v o l u m e ptn | Bcov pattern A pattern B−0.100.10.20.3 M ean c o rr e l a t i on i n v o l u m e ptn | BB pattern A pattern B−0.100.10.20.3 M ean c o rr e l a t i on i n v o l u m e ptn | Bcov+BB AR(1) A Figure 3. Analyses in real MRI data. Panels A-C: results obtained with an AR(1) for the lack of independence of sequential observations; panels D-F: same analyses obtained with FAST to model lack of independence. A and D: boxplots of mean volume correlations from the subjects in the sample. The boxplots on the left refer to RSA without partialling out confound terms; from second left to right, to RSA when partialling out
Bcov , BB , and both for patterns A and B (cf. Figure 1). B, E: t maps of RSA correlations without correction, showing large artefactual positive correlations for pattern A. C, F: searchlight analysis of the RSA correlation of the confounds and the data. Partialling out the off-diagonal terms of
Bcov removed a large part of global bias (panel A, ptn|Bcov , second from left), as did using BB ( ptn|BB , third from left). Indeed, Bcov and BB are similar (Figure 2). However, only using both Bcov and BB simultaneously achieved results comparable to those of the simulations ( ptn|Bcov+BB , right), suggesting the existence of local bias (i.e., bias that could not be captured by one covariance term alone). Testing this latter RSA model using a permutation method to correct for multiple testing on leave-one-out resampling of the original data gave a false positive rate of 0.07 for p < 0.05. To show where BB might have improved over Bcov in the combination, Panel C of Figure 3, right, displays the RSA correlation of BB itself after partialling out the effect of Bcov . The strongest correlations followed the distribution of vessel density (Viviani 2016), a main source of high variance in EPI data originating in blood pulsation (Lund et al. 2006). The increased signal on the midline is similar to effects of respiratory volumes reported by Birn et al. (2006). These voxels were missing from the RSA correlation of
Bcov , notwithstanding the high levels of signal noise. There was no difference when using
Svar instead of BB in these analyses. In panels D-F of Figure 3 we repeated the same analyses using FAST (Corbin et al. 2018), an option in the SPM package to model autocorrelations more flexibly by including several exponential time terms and their derivatives. The effectiveness of the correction was similar here as in the AR(1) model, suggesting the simpler model to be adequate in this case.
C F
Bcov BB|Bcov Bcov BB|Bcovr
B E pattern A pattern At -20 -10 0 -20 -10 0 -20 -10 0 -20 -10 0 -20 -10 0 -20 -10 0 pattern A pattern B−0.100.10.20.3 M ean c o rr e l a t i on i n v o l u m e ptn wo adjustment pattern A pattern B−0.100.10.20.3 ptn | Bcov pattern A pattern B−0.100.10.20.3 ptn | BB pattern A pattern B−0.100.10.20.3 ptn | Bcov+BB A pattern A pattern B−0.100.10.20.3 M ean c o rr e l a t i on i n v o l u m e ptn wo adjustment pattern A pattern B−0.100.10.20.3 ptn | Bcov pattern A pattern B−0.100.10.20.3 ptn | BB pattern A pattern B−0.100.10.20.3 ptn | Bcov+BB C FASTAR(1)
Figure 4. Analyses in real MRI data at short TR (same organization of panels as in Figure 3). Note that the scale of the t maps in B and E differs from the one used in Figure 2.
A question raised by these data is whether the effectiveness of
Bcov may be due by the ease with which the autocorrelation of residuals is modelled in these data. In Figure 4 we applied the same partial RSA strategy to correct for global bias in MRI data collected for the present study with TR = 0.34 sec (N = 32), where the correlation of sequential acquisitions may be larger. Perhaps surprisingly, the extent of global bias was smaller here than in the previous sample. Nevertheless, there were widespread artefactual RSA correlations (panels B and E), following the spatial pattern of
Bcov , which loaded above all on the center of the volume (where noise increases with distance from the coil) and to a lower extent on white matter (panel C). Importantly, the control for global bias through partial correlation was equally or more effective than in the long TR data (panels A and C). The use of FAST in this dataset (panel F) revealed the contribution of BB to be lower after adjusting for Bcov , as one would expect in a more realistic autocorrelation model of the error. There were no false positives in the leave-one-out resampling estimate of the FAST model after partialling out BB and Bcov at the p < 0.05 threshold, corrected at voxel level with the permutation method.
Additional strategies for the control of bias
Several additional strategies to control for bias are conceivable in individual studies. One is limiting the volume mask to gray matter to make the residual correlation model more representative of the target brain areas, or estimating this model in the pooled voxels where the RSA correlations are strongest. A more aggressive strategy consists in computing correlations between off-diagonal terms further away from the main diagonal, as those close to it are those where the effect of the design is strongest (Figure 2). This means omitting the cross-products of trials that are close to each other in the statistics, assuming that similarity can be estimated from further trials. In our dataset, implementing this strategy improved uncorrected global bias (Figure 5). However, partial correlation still appeared to improve the level of global bias in individual volumes. Figure 5. Estimates of global bias in the AR(1) model of long TR data when the sum of squares and cross products matrix is used with an off-diagonal off-set of 2 and 3 units.
Discussion
This work sought to address the confound arising in RSA from the non-orthogonality of the design matrix. Given the sluggishness of the BOLD response, it may not be easy to obtain orthogonality of regressors after convolving with the canonical BOLD response. Unfortunately, the precise form taken by the ensuing confound also depends on the temporal autocorrelation of the residuals in successive MR acquisitions. The importance of this issue in RSA contrasts with its benign consequences in ordinary analyses of functional imaging experiments. Although the imaging literature contains numerous statements on the consequences of inadequate models of residual temporal correlation on significance levels, starting at least from Purdon and Weisskoff (1998), these consequences only concern inference at the first level, which is almost never considered in practice. Irrespective of residual correlations and how they are modelled, coefficient estimates (including OLS estimates) are always unbiased (Diggle et al. 2002), and in ordinary functional imaging studies it is these coefficient maps that are brought to the second level, where they are independent. This fortunate state of affairs does not hold in RSA. There was no particular problem in our dataset to control for bias in data with short TR. The RSA correlation of of BB adjusted for Bcov shows that there was no effect here from large vessels, in contrast to the EPI data at longer TR. At this short TR, the movement covariates from realignment show oscillations at about 1Hz that are consistent with heartbeat. A plausible explanation of this finding is that the movement covariates added at the first level were capturing respiration and heartbeat variance, which may have generated high aliased autocorrelation in the data of the first sample (Lund et al. 2006). These results suggest that patterns of confound may depend on the distribution of noise when using different coils. Even when residual temporal autocorrelation is not adequately modelled, the ensuing global bias may be corrected by estimating the variance of the coefficients of the first-level model from the volume data. This was successful in our samples irrespective of the TR of the acquisition. It is more difficult to provide a comprehensive strategy to counter local bias. Our results suggest that local bias was minor in our data, as shown by the large effects of corrections obtained with single terms, but it would be difficult to generalize to all possible datasets. In the context of partial correlation, what one may do is to use several adjustment terms simultaneously, hoping to capture enough of the true variance-covariance varying across the volume. It seems a good idea to compare the effect of different correction terms and error correlation models, as true correlations with the similarity matrix of the stimuli should be robust to these variations. Different sets of correction terms may be applied to evaluate robustness of results. If there is variation in the patterns between individuals (as when the similarity matrices vary across the sample), it may be possible to use volume estimates of global bias (not centered) as a confounding covariate at second level, so that the coefficient of the constant term refers to a global bias of zero, and compare the outcome with the model without the confounding covariate. However, the possible existence of local bias implies no guarantee that any set of correction terms be equally and uniformly effective in all voxels. In this case, possible spurious associations with the similarity matrix of the stimuli may emerge. Methods
The study was based on the design of Labek et al. ( , which was chosen because it used a classic block design giving rise to strong bias. Participants viewed drawings of individuals, each displayed for 3 sec. in blocks of 4 drawings. Hence, each block had a duration of 12 sec. Scrambled drawings were displayed between blocks for the same duration to form a high-level baseline. There were three blocks of drawings depicting sad individuals and three blocks of control drawings. Drawings and blocks were assigned to arbitrary stimuli classes to form the patterns of Figure 1. Both EPI samples were acquired on the premises of the Clinic for Psychiatry and Psychotherapy of the University of Ulm, Germany, after obtaining written informed consent from volunteers recruited locally. The study was approved by the Ethical Committee of the University of Ulm. The simulations of Figure 1 were conducted in MATLAB. White noise volumes of isotropic voxel size 2mm were used as input data, creating 65 volumes for each of N = 30 subjects. In the analysis of Figure 2, data were collected with a T2*-weighted EPI in a Siemens 3T Allegra scanner (TR/TE 2260/35 msec, flip angle 90°, in-planar voxel size 3x3mm, slice thickness 3mm, 25% gap between slices, 35 slices). In each of N = 27 participants, 65 volumes were acquired. The data of Figure 4 were acquired in a Siemens 3T Prisma scanner with an EPI multiband sequence, with parameters chosen to minimize TR (slice acceleration factor 6, TR/TE 340/50 msec, flip angle 90°, in-planar voxel size 3x3mm, slice thickness 3mm, 13% gap between slices, 32 slices). No attempt was made to optimize for BOLD sensitivity, but we verified that the task and the contrast reported by Labek et al. ( were present in these data. In each of N = 32 participants, 433 images were acquired. To assist registration, a short T2-weighted EPI sequence of 5 volumes was acquired with no acceleration (TR/TE 4500/13 msec, flip angle 90°, same geometry as in the accelerated sequence). A 64-channels head coil was used with foam padding to minimize head motions. The EPI data were preprocessed with SPM12 (Functional Imaging Laboratory, London), using standard procedures (realignment, normalization), producing data resampled at the isotropic voxel size 2mm. Normalization of the multiband data was based on multimodal segmentation, which included the mean multiband and the mean EPI
The study was based on the design of Labek et al. ( , which was chosen because it used a classic block design giving rise to strong bias. Participants viewed drawings of individuals, each displayed for 3 sec. in blocks of 4 drawings. Hence, each block had a duration of 12 sec. Scrambled drawings were displayed between blocks for the same duration to form a high-level baseline. There were three blocks of drawings depicting sad individuals and three blocks of control drawings. Drawings and blocks were assigned to arbitrary stimuli classes to form the patterns of Figure 1. Both EPI samples were acquired on the premises of the Clinic for Psychiatry and Psychotherapy of the University of Ulm, Germany, after obtaining written informed consent from volunteers recruited locally. The study was approved by the Ethical Committee of the University of Ulm. The simulations of Figure 1 were conducted in MATLAB. White noise volumes of isotropic voxel size 2mm were used as input data, creating 65 volumes for each of N = 30 subjects. In the analysis of Figure 2, data were collected with a T2*-weighted EPI in a Siemens 3T Allegra scanner (TR/TE 2260/35 msec, flip angle 90°, in-planar voxel size 3x3mm, slice thickness 3mm, 25% gap between slices, 35 slices). In each of N = 27 participants, 65 volumes were acquired. The data of Figure 4 were acquired in a Siemens 3T Prisma scanner with an EPI multiband sequence, with parameters chosen to minimize TR (slice acceleration factor 6, TR/TE 340/50 msec, flip angle 90°, in-planar voxel size 3x3mm, slice thickness 3mm, 13% gap between slices, 32 slices). No attempt was made to optimize for BOLD sensitivity, but we verified that the task and the contrast reported by Labek et al. ( were present in these data. In each of N = 32 participants, 433 images were acquired. To assist registration, a short T2-weighted EPI sequence of 5 volumes was acquired with no acceleration (TR/TE 4500/13 msec, flip angle 90°, same geometry as in the accelerated sequence). A 64-channels head coil was used with foam padding to minimize head motions. The EPI data were preprocessed with SPM12 (Functional Imaging Laboratory, London), using standard procedures (realignment, normalization), producing data resampled at the isotropic voxel size 2mm. Normalization of the multiband data was based on multimodal segmentation, which included the mean multiband and the mean EPI without acceleration volumes simultaneously. The success of the segmentation was verified by inspection. Coefficients of noise and EPI data were modelled at the first level by convolving a stick function of duration 3 sec with a canonical BOLD curve. Each drawing presentation was modelled by its own regressor, giving as many beta images as stimuli (as suggested by Kriegeskorte et al. 2008). An intercept and realignment parameters (when applicable) were added as confounding covariates. No high-pass filter or autocorrelation model was applied to noise data. The high-pass filter of EPI data was set at 128 sec. No smoothing was applied to volume data prior to computing the first level model. The RSA was coded in MATLAB as an extension of the SPM package using the spm_searchlight function offered in that package. The spherical searchlight had a dimension of 8 mm, including 257 voxels when not located near the edge of the selected volume, defined by a mask including grey and white matter as well as CSF using the prior provided by SPM for these compartments at threshold p < 0.5. Searchlight volumes including less than 27 voxels were excluded. The similarity matrix contained the value zero for entries of different categories, and a constant value for entries of the same categories. As described in the main text, the similarity of the brain response was encoded as the sum of squares and cross products of the coefficients of the design shown in Figure 1. Pearson correlations and partial correlations between the upper diagonal terms of these two matrices were computed using the native MATLAB functions. Bcov was extracted from the internal SPM computation of that quantity, whereas
Scov and BB were computed as indicated by the equations in the main text. The overlays of panels B and E of Figures 3 and 4 were obtained by computing t maps of one-sample t test of the searchlight correlations brought to the second level and smoothed with a Gaussian kernel (FWHM 4mm). The overlays of panels C and F are average maps of searchlight correlations (same smoothing). The false positive rates reported in text were computed by resampling the data without replacement, each time leaving out one individual at the second level, and computing significance levels with a permutation method for strong control of significance levels (percentiles of maximal t values over the volume, 2000 permutations per resampling; see Holmes et al. 1996 ). The rates are the percent of resamples where the null was rejected at the voxel level-corrected significance threshold p < 0.05.
Acknowledgments
This work was conducted within the framework of the “Austrian NeuroCloud”, supported by the Austrian Federal Ministry of Education, Science and Research. Acquisition of data was supported by a collaborative grant from the Federal Institute for Drugs and Medical Devices (BfArM, Bonn, Grant No. V-17568/68502/2017-2020).
References
Birn, R.M., Diamond, J.B., Smith, M.A., Bandettini, P.A., 2006. Separating respiratory-variation-related fluctuations from neuronal-activity-related fluctuations in fMRI. NeuroImage 31, 1536-1548. Cai, M.B., Schuck, N.W., Pillow, J.W., Niv, Y., 2019. Representational structure of task structure? Bias in neural representational similarity analysis and a Byesian method for reducing bias. PLoS Comp. Biol. 15(5), e1006299.