（Group) Bayesian Representational Similarity Analysis

Overview

This demo shows how to use the Bayesian Representational Similarity Analysis (Bayesian RSA, Cai et al., 2016,2019) with a simulated dataset.

Bayesian RSA builds a generative model for fMRI data which include a covariance structure as estimation target. The covariance structure specifies the distribution by which the unknown spatial neural response patterns for each task condition follow. In the model, the task-evoked patterns modulated by the design matrix, together with spontaneous activity and noise, constitute the fMRI data. By marginalizing the intermediate unknown variables including voxel-wise response amplitudes to the task conditions, the model computes a log likelihood for any possible covariance structure of condition-evoked response patterns to give rise to the fMRI data. It can then directly infer the most likely covariance structure, and convert it to the similarity structure of the neural representation for the task conditions in a region of interest. This is different from the traditional RSA approach which first estimates neural response patterns from data and then treats the noisy estimated patterns as true patterns to calculate their similarity structure. It was shown (Cai et al., 2016,2019) that the traditional approach can introduce bias to the estimated similarity structure and BRSA reduces this bias.

The BrainIAK Bayesian RSA module (brainiak.reprsimil.brsa) includes the following variants of Bayesian RSA:

• (the original) BRSA: a model that estimates representational similarity structure for one subject, and maximum likelihood estimates of SNR and noise parameters for each voxel.
• GBRSA: a model that can either estimate representaitonal similarity structure for one subject, or the shared similarity structure across a group of subjects. It marginalizes voxel-wise SNR and noise parameters and is likely more accurate, but may be slower than BRSA. "G" stands for "Group".

This notebook demonstrates the usage of GBRSA on multiple subjects using simulated data. The usage of BRSA method is similar, and readers can refer to this example

To run GBRSA on single subject, you can provide data to GBRSA in the same format as for a group of subjects shown in this notebook, except that the data is a list of numpy array of one subject's data

Annotated Bibliography

1. Cai, M. B., Schuck, N. W., Pillow, J. W., & Niv, Y. (2016). A Bayesian method for reducing bias in neural representational similarity analysis. Advances in Neural Information Processing Systems (pp. 4951-4959). link Describes potential biases in computing RSA that could lead to spurious results and introduces a Bayesian approach to reduce bias.

2. Cai, M. B., Schuck, N. W., Pillow, J. W., & Niv, Y. (2019). Representational structure or task structure? Bias in neural representational similarity analysis and a Bayesian method for reducing bias. PLoS Computational Biology, 15(5), e1006299. link This paper improves the earlier version of BRSA by additionally modeling the spatial noise correlation and marginalizing voxel-wise noise parameters. The paper also introduces group BRSA and the use of cross-validation to compare an estimated model against a null model, and further extends BRSA to task-signal decoding, using the estimated similarity structure as an empirical prior for estimating neural patterns.

Table of contents

• Simulate data
• Fit Group Bayesian RSA to simulated data
• Evaluation
• Other usage: "decoding" task-related signal from new data
• Model selection by cross-validataion
• Avoid false discovery
• Summary

Notes on data preparation

When you apply this tool to real fMRI data, it is required that the data of each participant to be motion corrected. If multiple runs are acquired for each participant, they should be spatially aligned. Slice-timing correction is recommended.

You will need to have the mask of the Region of Interest (ROI) ready (typically defined anatomically or by independent tasks). nilearn provides tools to extract signal from mask. You can refer to this tutorial

Notes on model assumption

Please note that the model assumes that the covariance matrix U which all $\beta_i$ follow describe a multi-variate Gaussian distribution that is zero-meaned. This assumption does not imply that there must be both positive and negative responses across voxels. However, it means that (Group) Bayesian RSA treats the task-evoked activity against baseline BOLD level as signal, while in other RSA tools the deviation of task-evoked activity in each voxel from the average task-evoked activity level across voxels may be considered as signal of interest. Due to this assumption in (G)BRSA, relatively high degree of similarity may be expected when the activity patterns of two task conditions share a strong sensory driven components. When two task conditions elicit exactly the same activity pattern but only differ in their global magnitudes, under the assumption in (G)BRSA, their similarity is 1. However, if one take the assumption that only deviation of pattern from average patterns is signal of interest (which is currently not supported by (G)BRSA), their similarity would be -1 because the deviations of the two patterns from their average pattern are exactly opposite.

Load some package which we will use in this demo.

If you see error related to loading any package, you can install that package. For example, if you use Anaconda, you can use "conda install matplotlib" to install matplotlib.

%matplotlib inline
import scipy.spatial.distance as spdist
from scipy import stats
import numpy as np
from brainiak.reprsimil.brsa import GBRSA
import brainiak.utils.utils as utils
import matplotlib.pyplot as plt
import logging
np.random.seed(10)

You might want to keep a log of the output.

logging.basicConfig(
    level=logging.DEBUG,
    filename='gbrsa_example.log',
    format='%(relativeCreated)6d %(threadName)s %(message)s')

Simulate data

We want to simulate some data in which each voxel responds to different task conditions differently, but following a common covariance structure.

This example simulate 5 subjects. If you find the whole notebook runs for too long on your computer, you can reduce the number of simulated subjects by changing n_subj below.

Load an example design matrix.

To use BRSA, you need to prepare design matrix for the task with your favorate software, such as using 3ddeconvolve of AFNI, or using SPM or FSL.

The design matrix reflects your belief of how fMRI signal should respond to a task (if a voxel does respond). The common assumption is that a neural event that you are interested in will elicit a slow hemodynamic response in some voxels. The response peaks around 4-6 seconds after the event onset and dissipatess ~12 seconds after the event. Therefore, typically you convolve a time series A, composed of delta (stem) functions reflecting the time of each neural event belonging to the same category (e.g. all trials in which a participant sees a face), with a hemodynamic response function B, to form the hypothetic response of any voxel to such type of neural event. One convoluted time course can be generated this way for each type of event. The time courses corresponding to all events, put together, are called design matrix. Our goal is to figure out how the (spatial) response patterns of a population of voxels (in an Region of Interest, ROI) to different types of tasks (e.g., watching different categories of animals, different conditions of a cognitive task) are similar or disimilar. So we need the design matrix in order to estimate the similarity matrix we are interested in.

Basically, (G)BRSA just needs a numpy array which is in size of {time points} * {condition}. We can use the utility called ReadDesign from brainiak.utils to read a design matrix generated from AFNI. For design matrix saved as Matlab data file by SPM or or other toolbox, you can use scipy.io.loadmat('YOURFILENAME') and extract the design matrix from the dictionary returned. You can also generate design matrix using the function gen_design in brainiak.utils. It takes in (names of) event timing files in AFNI or FSL format (denoting onsets, duration, and weight for each event belonging to the same condition) and outputs the design matrix as numpy array.

details about design matrix

In this simulation, we use a design matrix of a task consisting of 16 different conditions, which by design cannot be fully counter-balanced in their orders. For simplicity, we repeat the design matrix of one run by 2 to 3 times, mimicking 2 to 3 fMRI runs with identical timing. Note that different subjects do not have to have the same number of voxels or time points. The timing of the task conditions can also differ across subjects. The simulation below reflects this.

In typical fMRI analysis, some nuisance regressors such as head motion, baseline time series and slow drift are also entered into regression. In using (G)BRSA, you should not include such nuisance regressors into the design matrix, because the spatial spread of such nuisance regressors might be quite different from the spatial spread of task related signal. Including such nuisance regressors in design matrix might influence the pseudo-SNR map, which in turn influences the estimation of the shared covariance matrix. But you may include motion time course in the nuisance parameter.

n_subj = 5
n_run = np.random.randint(2, 4, n_subj)
ROI_edge = np.random.randint(20, 40, n_subj)
# We simulate "ROIs" of a square shape
design = [None] * n_subj
for subj in range(n_subj):
    design[subj] = utils.ReadDesign(fname="example_design.1D")
    design[subj].n_TR = design[subj].n_TR * n_run[subj]
    design[subj].design_task = np.tile(design[subj].design_task[:,:-1],
                                 [n_run[subj], 1])
    # The last "condition" in design matrix
    # codes for trials subjects made an error.
    # We ignore it here.
n_C = np.size(design[0].design_task, axis=1)
# The total number of conditions.
n_V = [int(roi_e**2) for roi_e in ROI_edge]
# The total number of simulated voxels
n_T = [d.n_TR for d in design]
# The total number of time points,
# after concatenating all fMRI runs
fig = plt.figure(num=None, figsize=(12, 3),
                 dpi=150, facecolor='w', edgecolor='k')
plt.plot(design[0].design_task)
plt.ylim([-0.2, 0.4])
plt.title('hypothetic fMRI response time courses '
          'of all conditions for one subject\n'
         '(design matrix)')
plt.xlabel('time')
plt.show()

simulate noise

We simulate noise which is Gaussian Process (mimicking smooth noise in fMRI) in space and AR(1) in time

noise_bot = 0.5
noise_top = 1.5
noise_level = [None] * n_subj
noise = [None] * n_subj
rho1 = [None] * n_subj
white_noise_magnitude = 0.4
for subj in range(n_subj):
    # each subject may have different noise level
    noise_level[subj] = np.random.rand(n_V[subj]) * \
        (noise_top - noise_bot) + noise_bot
    # The standard deviation of the noise is in the range of [noise_bot, noise_top]
    # In fact, we simulate autocorrelated noise with AR(1) model. So the noise_level reflects
    # the independent additive noise at each time point (the "fresh" noise)

# AR(1) coefficient
rho1_top = 0.8
rho1_bot = -0.1
for subj in range(n_subj):
    rho1[subj] = np.random.rand(n_V[subj]) \
        * (rho1_top - rho1_bot) + rho1_bot

noise_smooth_width = 10.0
dist2 = [None] * n_subj
for subj in range(n_subj):
    coords = np.mgrid[0:ROI_edge[subj], 0:ROI_edge[subj], 0:1]
    coords_flat = np.reshape(coords,[3, n_V[subj]]).T
    dist2[subj] = spdist.squareform(spdist.pdist(coords_flat, 'sqeuclidean'))

    # generating noise
    K_noise = noise_level[subj][:, np.newaxis] \
        * (np.exp(-dist2[subj] / noise_smooth_width**2 / 2.0) \
           + np.eye(n_V[subj]) * white_noise_magnitude ** 2) * noise_level[subj]
    # We make spatially correlated noise by generating
    # noise at each time point from a Gaussian Process
    # defined over the coordinates.
    L_noise = np.linalg.cholesky(K_noise)
    noise[subj] = np.zeros([n_T[subj], n_V[subj]])
    noise[subj][0, :] = np.dot(L_noise, np.random.randn(n_V[subj]))\
        / np.sqrt(1 - rho1[subj]**2)
    for i_t in range(1, n_T[subj]):
        noise[subj][i_t, :] = noise[subj][i_t - 1, :] * rho1[subj] \
            + np.dot(L_noise,np.random.randn(n_V[subj]))
    # For each voxel, the noise follows AR(1) process:
    # fresh noise plus a dampened version of noise at
    # the previous time point.
    # In this simulation, we also introduced spatial smoothness resembling a Gaussian Process.
    # Notice that we simulated in this way only to introduce spatial noise correlation.
    # Gaussian Process is not the assumption of the form of spatial noise correlation in the model.
    noise[subj] += np.random.randn()
plt.pcolor(K_noise)
plt.colorbar()
plt.xlim([0, ROI_edge[-1] * ROI_edge[-1]])
plt.ylim([0, ROI_edge[-1] * ROI_edge[-1]])
plt.title('Spatial covariance matrix of noise\n of one participant in a simulated square "ROI"')
plt.xlabel('x coordinates of voxels')
plt.ylabel('y coordinates of voxels')
plt.show()

fig = plt.figure(num=None, figsize=(12, 2), dpi=150,
                 facecolor='w', edgecolor='k')
plt.plot(noise[-1][:, 0])
plt.title('noise in an example voxel')
plt.xlabel('TR')
plt.ylabel('demeaned signal')
plt.show()