brainiak.reprsimil package

A Bayesian method to perform Representational Similarity Analysis


brainiak.reprsimil.brsa module

Bayesian Representational Similarity Analysis (BRSA)

This implementation is based the work in [Cai2016].
[Cai2016]“A Bayesian method for reducing bias in neural representational similarity analysis”, M.B. Cai, N. Schuck, J. Pillow, Y. Niv, Neural Information Processing Systems (NIPS) 29, 2016.
class brainiak.reprsimil.brsa.BRSA(n_iter=50, rank=None, GP_space=False, GP_inten=False, tol=0.002, auto_nuisance=True, n_nureg=6, verbose=False, eta=0.0001, space_smooth_range=None, inten_smooth_range=None, tau_range=10.0, init_iter=20, optimizer=’BFGS’, rand_seed=0)

Bases: sklearn.base.BaseEstimator

Bayesian representational Similarity Analysis (BRSA)

Given the time series of neural imaging data in a region of interest (ROI) and the hypothetical neural response (design matrix) to each experimental condition of interest, calculate the shared covariance matrix of the voxels(recording unit)’ response to each condition, and the relative SNR of each voxels. The relative SNR could be considered as the degree of contribution of each voxel to this shared covariance matrix. A correlation matrix converted from the covariance matrix will be provided as a quantification of neural representational similarity.

\[ \begin{align}\begin{aligned}Y = X \cdot \beta + \epsilon\\\beta_i \sim N(0,(s_{i} \sigma_{i})^2 U)\end{aligned}\end{align} \]
  • n_iter (int, default: 200) – Number of maximum iterations to run the algorithm.
  • rank (int, default: None) – The rank of the covariance matrix. If not provided, the covariance matrix will be assumed to be full rank. When you have many conditions (e.g., calculating the similarity matrix of responses to each event), you might want to start with specifying a lower rank and use metrics such as AIC or BIC to decide the optimal rank.
  • auto_nuisance (boolean, default: True) – In order to model spatial correlation between voxels that cannot be accounted for by common response captured in the design matrix, we assume that a set of time courses not related to the task conditions are shared across voxels with unknown amplitudes. One approach is for users to provide time series which they consider as nuisance but exist in the noise (such as head motion). The other way is to take the first n_nureg principal components in the residual after subtracting the response to the design matrix from the data, and use these components as the nuisance regressor. If this flag is turned on, the nuisance regressor provided by the user is used only in the first round of fitting. The PCs from residuals will be used in the next round of fitting. Note that nuisance regressor is not required from user. If it is not provided, DC components for each run will be used as nuisance regressor in the initial fitting.
  • n_nureg (int, default: 6) – Number of nuisance regressors to use in order to model signals shared across voxels not captured by the design matrix. This parameter will not be effective in the first round of fitting.
  • GP_space (boolean, default: False) – Whether to impose a Gaussion Process (GP) prior on the log(pseudo-SNR). If true, the GP has a kernel defined over spatial coordinate. This is experimental and slow. I find that when SNR is genrally low, smoothness can be overestimated. But I think this is better than not imposing any regularization.
  • GP_inten (boolean, defualt: False) – Whether to include a kernel definved over the intensity of image. GP_space should be True as well if you want to use this, because the smoothness should be primarily in space. Smoothness in intensity is just complementary.
  • tol (tolerance parameter passed to the minimizer.) –
  • verbose (boolean, default: False) – Verbose mode flag.
  • eta (a small number added to the diagonal element of the) – covariance matrix in the Gaussian Process prior. This is to ensure that the matrix is invertible.
  • space_smooth_range (the distance (in unit the same as what) – you would use when supplying the spatial coordiates of each voxel, typically millimeter) which you believe is the maximum range of the length scale parameter of Gaussian Process defined over voxel location. This is used to impose a half-Cauchy prior on the length scale. If not provided, the program will set it to half of the maximum distance between all voxels.
  • inten_smooth_range (the difference in image intensity which) – you believe is the maximum range of plausible length scale for the Gaussian Process defined over image intensity. Length scales larger than this are allowed, but will be penalized. If not supplied, this parameter will be set to half of the maximal intensity difference.
  • tau_range (the reasonable range of the standard deviation) – of the Gaussian Process. Since the Gaussian Process is imposed on the log(SNR), this range should not be too large. 10 is a pretty loose range. This parameter is used in a half-Cauchy prior on the standard deviation
  • init_iter (how many initial iterations to fit the model) – without introducing the GP prior before fitting with it, if GP_space or GP_inten is requested. This initial fitting is to give the parameters a good starting point.
  • optimizer (the optimizer to use for minimizing cost function.) – We use ‘BFGS’ as a default. Users can try other optimizer coming with scipy.optimize.minimize, or a custom optimizer.
  • rand_seed (int, default: 0) – Seed for initializing the random number generator.

The shared covariance matrix, shape=[condition,condition].


The Cholesky factor of the shared covariance matrix – (lower-triangular matrix), shape=[condition,condition].


the correlation matrix derived from the shared covariance matrix, – shape=[condition,condition]


array, shape=[voxels,] – The pseuso-SNR of all voxels. They are normalized such that the geometric mean is 1


array, shape=[voxels,] – The estimated standard deviation of the noise in each voxel Assuming AR(1) model, this means the standard deviation of the refreshing noise.


array, shape=[voxels,] – The estimated autoregressive coefficient of each voxel


scalar, only if GP_space or GP_inten is True. – the standard deviation of the GP prior


scalar, only if GP_space or GP_inten is True – the length scale of Gaussian Process prior of log(SNR)


scalar, only if GP_inten is True – the length scale in fMRI intensity of the GP prior of log(SNR)


array, shape=[conditions, voxels] – The maximum a posterior estimation of the response amplitudes of each voxel to each task condition.


array, shape=[n_nureg, voxels] – The loading weights of each voxel for the shared time courses not captured by the design matrix.

fit(X, design, nuisance=None, scan_onsets=None, coords=None, inten=None)

Compute the Bayesian RSA

  • X (2-D numpy array, shape=[time_points, voxels]) – If you have multiple scans of the same participants that you want to analyze together, you should concatenate them along the time dimension after proper preprocessing (e.g. spatial alignment), and specify the onsets of each scan in scan_onsets.
  • design (2-D numpy array, shape=[time_points, conditions]) – This is the design matrix. It should only include the hypothetic response for task conditions. You do not need to include regressors for a DC component or motion parameters, unless with a strong reason. If you want to model head motion, you should include them in nuisance regressors. If you have multiple run, the design matrix of all runs should be concatenated along the time dimension, with one column across runs for each condition.
  • nuisance (optional, 2-D numpy array,) – shape=[time_points, nuisance_factors] The responses to these regressors will be marginalized out from each voxel, which means they are considered, but won’t be assumed to share the same pseudo-SNR map with the design matrix. Therefore, the pseudo-SNR map will only reflect the relative contribution of design matrix to each voxel. You can provide time courses such as those for head motion to this parameter. Note that if auto_nuisance is set to True, this input will only be used in the first round of fitting. The first n_nureg principal components of residual (excluding the response to the design matrix) will be used as the nuisance regressor for the second round of fitting. If auto_nuisance is set to False, the nuisance regressors supplied by the users together with DC components will be used as nuisance time series.
  • scan_onsets (optional, an 1-D numpy array, shape=[runs,]) – This specifies the indices of X which correspond to the onset of each scanning run. For example, if you have two experimental runs of the same subject, each with 100 TRs, then scan_onsets should be [0,100]. If you do not provide the argument, the program will assume all data are from the same run. The effect of them is to make the inverse matrix of the temporal covariance matrix of noise block-diagonal.
  • coords (optional, 2-D numpy array, shape=[voxels,3]) – This is the coordinate of each voxel, used for implementing Gaussian Process prior.
  • inten (optional, 1-D numpy array, shape=[voxel,]) – This is the average fMRI intensity in each voxel. It should be calculated from your data without any preprocessing such as z-scoring. Because it should reflect whether a voxel is bright (grey matter) or dark (white matter). A Gaussian Process kernel defined on both coordinate and intensity imposes a smoothness prior on adjcent voxels but with the same tissue type. The Gaussian Process is experimental and has shown good performance on some visual datasets.