Navigating fMRI analysis techniques: a practical guide

Welcome to Navigating fMRI analysis techniques: a practical guide!

How would you like to begin?

  • Interactive flowchart (in-progress)
  • Non-interactive flowchart (pdf)
  • Full list of fMRI analysis techniques
  • About this tool
  • Contact: scottibrain@gmail.com

    PLEASE NOTE THAT THIS WEBSITE IS NOT FINISHED. PLEASE CONTACT AUTHORS FOR INFORMATION.

    NavigateFMRI (or, Navigating fMRI analysis techniques: a practical guide) was developed by Paul S. Scotti, Jiageng Chen, Xiaoli Zhang, and Julie D. Golomb to help researchers to better understand the currently available fMRI analysis techniques for exploring fMRI brain data and testing hypotheses about the brain.

    More specifically, the goal of this tool is to review the most widely used post-processing fMRI methods, including both basic and advanced computational techniques, in an intuitive and hands-on manner. We provide a flowchart that guides the reader to the most appropriate analysis methods, and then we supply simulated data and example Python code to implement each method. Unlike other fMRI resources (we highly recommend Brainiak and fmri4newbies), this web tool focuses on the big picture--what are the different methods, when would you use each method, how do the methods relate to each other and the underlying research question, and what are the basic steps involved in each method (i.e., breadth over depth).

    There are many stages to running a complete fMRI experiment, including experimental design, data collection, data pre-processing (e.g., motion correction), data processing (e.g., general linear model estimation), and finally post-processing analysis (e.g., univariate analysis) to reach a conclusion about the brain. Due to the number of steps prior to post-processing analysis, the teaching of these methods can often become confusing despite relatively straightforward steps. For instance, functional connectivity is just the correlations of time series, and representational similarity analysis is just the correlations of conditions! Such intuitive understandings can often be lost because when a researcher finally reaches the post-processing stage, they may not have easy access to their data content (e.g., interpreting the contents of an SPM.mat file is a cruel fate) and only be able to access it via package- and method-specific toolboxes that further obscure our understanding of the basic steps performed under the hood.

    This tool can help readers to better understand the main ideas and core implementation of a variety of basic and advanced fMRI methods without worrying about the other fMRI stages such as pre-processing. We accomplish this through easy-to-follow flowcharts, high-level summaries of each method, and practical Python examples using simulated data (we only use standard Python packages in our examples--no custom datasets, packages, or toolboxes need to be installed).

    Go back

    What kind of research question are you trying to answer?

    Work in progress

    Go back

    List of fMRI analysis techniques

    note: all below pages are works in progress!
  • Univariate analysis
  • Functional connectivity
  • Searchlight maps / T- and F-maps
  • Repetition suppression
  • Multivoxel pattern analysis: correlation
  • Multivoxel pattern classifier
  • Representational similarity analysis
  • Dimensionality reduction
  • Encoding & decoding model
  • Phase-encoded mapping / Traveling wave
  • Population receptive field mapping
  • Hidden markov model
  • Neural networks
  • Graph theory
  • Effective connectivity / Dynamic causal model
  • Mixed effects modeling
  • Go back

    Univariate Analysis

    For a univariate analysis, you need, for each participant, the brain activations for one condition and the average brain activations for another condition (raw BOLD signal or beta weights can be used as the estimate of brain activity, where beta weights are obtained from general linear model estimation) within a brain region. Then you can simply conduct a t-test across participants on the difference in brain activations between conditions. This allows you to test whether a brain region is consistently more or less active across people in response to certain conditions.

    Upsides

    lorem

    Downsides

    lorem

    Example

    First import the data and a few standard Python packages.
    
            # general packages
            import numpy as np
            import matplotlib.pyplot as plt
            import scipy as sp
            
            # pull data from GitHub
            import requests, io
            for array in ['activations','sub_id','cond_id','vox_id']:
                globals()['{}'.format(array)] = np.load(io.BytesIO(requests.get(
                    'http://paulscotti.github.io/NavigateFMRI/methods/univariate/{}.npy'.format(array)).content))
            
            # overview of the data
            num_subjects = len(np.unique(sub_id)) #20
            num_conditions = len(np.unique(cond_id)) #2
            num_voxels = len(np.unique(vox_id)) #50
            
    We have 20 subjects each with a 50-voxel brain region of interest (ROI). We have activations for every voxel for each of two conditions (horror movie and comedy movie).
    
            plt.imshow(np.reshape(activations[cond_id==0],(num_subjects,num_voxels)),cmap='PuBu')
            plt.colorbar()
            plt.title("Activations for comedy movies")
            plt.yticks([]);plt.ylabel('Subject ID')
            plt.xticks([]);plt.xlabel('Voxel ID')
            plt.show()
            
            plt.imshow(np.reshape(activations[cond_id==1],(num_subjects,num_voxels)),cmap='OrRd')
            plt.colorbar()
            plt.title("Activations for horror movies")
            plt.yticks([]);plt.ylabel('Subject ID')
            plt.xticks([]);plt.xlabel('Voxel ID')
            plt.show()
            

    Let's average the activations across voxels for each subject, for each condition. Each cell represents a single subject's average activation. Then take the differences in activations between conditions for each subject.
    
            plt.imshow(np.array([np.mean(np.reshape(activations[cond_id==0],(num_subjects,num_voxels)),axis=1)]).T,cmap='PuBu',clim=[0,1])
            plt.title("Avg")
            plt.yticks([]);plt.xticks([]);plt.colorbar()
            plt.show()
        
            plt.imshow(np.array([np.mean(np.reshape(activations[cond_id==1],(num_subjects,num_voxels)),axis=1)]).T,cmap='OrRd',clim=[0,1])
            plt.title("Avg")
            plt.yticks([]);plt.xticks([]);plt.colorbar()
            plt.show()
        
            diffs=np.array([np.mean(np.reshape(activations[cond_id==1],(num_subjects,num_voxels)),axis=1)]).T - np.array([np.mean(np.reshape(activations[cond_id==0],(num_subjects,num_voxels)),axis=1)]).T
            plt.imshow(diffs,cmap='bwr',clim=[-1,1])
            plt.title("Diff.")
            plt.yticks([]);plt.xticks([])
            plt.colorbar()
            plt.show()
            

    We can see that perhaps activations appear greater for comedy (blue) compared to red (horror). But is this statistically significant? To test this, let us now conduct a paired t-test between this array of differences and an array of zeros (if p-value is less than .05, can conclude significance). (Note: you could skip a step and simply do a paired t-test between the average activations per condition).
    
                stat = sp.stats.ttest_rel(diffs.flatten(),np.zeros(num_subjects))
                print("Horror vs. Comedy (paired t-test):\nt={:.4f}, p={:.4f}".format(stat.statistic,stat.pvalue))
            
    Horror vs. Comedy (paired t-test): t=1.7043, p=0.1046

    Relevant resources / publications

    Go back

    Repetition suppression

    Repetition suppression, or fMRI adaptation, refers to how neural activity adapts to the repeated presentation of the same stimulus. For example, if you present an image of Jennifer Aniston, the fusiform face area (FFA) will show increased activation. If you then immediately present the same image again, activation will not be as pronounced as the first presentation because the FFA specializes in the processing of faces and would adapt to the repeated face image. Alternatively, if you showed a picture of Barack Obama, activation may be as high as the first presentation of Jennifer Aniston because the FFA had adapted to Jennifer Aniston and not Barack Obama. This difference can then be contrasted with a region of the brain that is not selective for faces, such as early visual cortex, where such adaptation would likely not be observed. By testing for reduced activation to repeated presentations of a stimulus, one can gauge the selectivity of a brain region.

    Upsides

    lorem

    Downsides

    lorem

    Example

    Please see the univariate page for an example. The difference here would be that the conditions would consist of same-genre vs. different-genre trial repetitions.

    Relevant resources / publications

  • Dilks, D. D., Julian, J. B., Kubilius, J., Spelke, E. S., & Kanwisher, N. (2011). Mirror-image sensitivity and invariance in object and scene processing pathways. Journal of Neuroscience, 31(31), 11305-11312.
  • Go back

    Searchlight / T-map / F-map

    If the goal of the researcher is to explore where the brain is most responsive to a task, it may be useful to create a statistical map (t-map, f-map) or to employ searchlight mapping. Essentially, a statistical map involves repeating the analysis for every voxel separately and plotting the corresponding statistical value (e.g., t or F value) for each voxel onto the brain. A searchlight map involves repeated analyses for a number of spherical regions of interest. Imagine dividing up your entire brain into many spherical regions and then likewise calculating the corresponding t-value for every sphere.

    Upsides

    lorem

    Downsides

    lorem

    Example of searchlight mapping

    First import the data and a few standard Python packages.
    
            # general packages
            import numpy as np
            import matplotlib.pyplot as plt
            import scipy as sp
            
            # pull data from GitHub
            import requests, io
            for array in ['activations1','activations2','coord']:
                globals()['{}'.format(array)] = np.load(io.BytesIO(requests.get(
                    'http://paulscotti.github.io/NavigateFMRI/methods/searchlight/{}.npy'.format(array)).content))
            # activations1 are activations for comedy movies
            # activations2 are activations for horror movies
            
    Now let's visualize our brain, composed of voxels with each voxel representing a certain activation strength (e.g., beta weights from general linear model). Note that this simulated brain is a cube for simplicity, where the x y and z axes all span a range of 20 voxels. This means that there are 8,000 total voxels in this brain-cube (20 cubed = 8,000). We have voxel activations for each of two conditions: one where the participant was watching a comedy movie and another where the participant was watching a horror movie.
    
            from mpl_toolkits.mplot3d import Axes3D
            fig = plt.figure()
            ax = fig.gca(projection='3d')
            p = ax.scatter(coord[:,0],
                        coord[:,1],
                        coord[:,2],
                        c=activations1,
                        marker='s',vmin=0,vmax=2,
                        cmap='PuBu')
            fig.colorbar(p)
            plt.title("Activations for comedy movies")
            plt.show()
            
            fig = plt.figure()
            ax = fig.gca(projection='3d')
            p = ax.scatter(coord[:,0],
                        coord[:,1],
                        coord[:,2],
                        c=activations1,
                        marker='s',vmin=0,vmax=2,
                        cmap='OrRd')
            fig.colorbar(p)
            plt.title("Activations for horror movies")
            plt.show()
            

    Now let's pick a random location and cut out a sphere with a radius of 5 voxels for each of the above brain-cubes, and calculate the difference.
    
            x0 = 12
            y0 = 12
            z0 = 12
            center=np.array([x0,y0,z0])
            cutoff=5
            
            distance=sp.spatial.distance.cdist(coord,center.reshape(1,-1)).ravel()
            points_in_sphere=coord[distance<cutoff]
            
            fig = plt.figure()
            ax = fig.gca(projection='3d')
            p = ax.scatter(points_in_sphere[:,0],
                        points_in_sphere[:,1],
                        points_in_sphere[:,2],
                        c=activations1[distance<cutoff],
                        marker='s',cmap='PuBu',vmin=0,vmax=2)
            fig.colorbar(p)
            ax.set_xlim3d(0,20)
            ax.set_ylim3d(0,20)
            ax.set_zlim3d(0,20)
            plt.title("Comedy")
            plt.show()
            
            fig = plt.figure()
            ax = fig.gca(projection='3d')
            p = ax.scatter(points_in_sphere[:,0],
                        points_in_sphere[:,1],
                        points_in_sphere[:,2],
                        c=activations2[distance<cutoff],
                        marker='s',cmap='OrRd',vmin=0,vmax=2)
            fig.colorbar(p)
            ax.set_xlim3d(0,20)
            ax.set_ylim3d(0,20)
            ax.set_zlim3d(0,20)
            plt.title("Horror")
            plt.show()
            
            fig = plt.figure()
            ax = fig.gca(projection='3d')
            p = ax.scatter(points_in_sphere[:,0],
                        points_in_sphere[:,1],
                        points_in_sphere[:,2],
                        c=activations2[distance<cutoff]-activations1[distance<cutoff],
                        marker='s',cmap='bwr',vmin=-3,vmax=3)
            fig.colorbar(p)
            ax.set_xlim3d(0,20)
            ax.set_ylim3d(0,20)
            ax.set_zlim3d(0,20)
            plt.title("Difference")
            plt.show()
            

    Let's repeat this procedure for every possible sphere location, keeping track of the average difference between conditions for every sphere.
    
            diffs=np.full(num_voxels,np.nan)
            coord_tracker=np.full((num_voxels,3),np.nan)
            cnt=0
            for x0 in range(20):
                for y0 in range(20):
                    for z0 in range(20):
                        center=np.array([x0,y0,z0])
                        distance=sp.spatial.distance.cdist(coord,center.reshape(1,-1)).ravel()
                        points_in_sphere=coord[distance<cutoff]
                        diffs[cnt] = np.abs(np.mean(activations2[distance<cutoff]-activations1[distance<cutoff]))
                        coord_tracker[cnt,:]=[x0,y0,z0]
                        cnt+=1
    
            print("Largest absolute difference found at x={} y={} z={}".format(coord_tracker[np.argmax(diffs),0],
                                                                        coord_tracker[np.argmax(diffs),1],
                                                                        coord_tracker[np.argmax(diffs),2]))
            

    Largest absolute difference found at x=6.0 y=6.0 z=6.0

    And now let's take a look at the sphere where there was the largest difference of conditions! In this way we have explored where the brain was most responsive to task-specific information.
    
            x0 = 6
            y0 = 6
            z0 = 6
            center=np.array([x0,y0,z0])
    
            distance=sp.spatial.distance.cdist(coord,center.reshape(1,-1)).ravel()
            points_in_sphere=coord[distance<cutoff]
    
            fig = plt.figure()
            ax = fig.gca(projection='3d')
            p = ax.scatter(points_in_sphere[:,0],
                    points_in_sphere[:,1],
                    points_in_sphere[:,2],
                    c=activations1[distance<cutoff],
                    marker='s',cmap='PuBu',vmin=0,vmax=2)
            fig.colorbar(p)
            ax.set_xlim3d(0,20)
            ax.set_ylim3d(0,20)
            ax.set_zlim3d(0,20)
            plt.title("Comedy")
            plt.show()
    
            fig = plt.figure()
            ax = fig.gca(projection='3d')
            p = ax.scatter(points_in_sphere[:,0],
                    points_in_sphere[:,1],
                    points_in_sphere[:,2],
                    c=activations2[distance<cutoff],
                    marker='s',cmap='OrRd',vmin=0,vmax=2)
            fig.colorbar(p)
            ax.set_xlim3d(0,20)
            ax.set_ylim3d(0,20)
            ax.set_zlim3d(0,20)
            plt.title("Horror")
            plt.show()
    
            fig = plt.figure()
            ax = fig.gca(projection='3d')
            p = ax.scatter(points_in_sphere[:,0],
                    points_in_sphere[:,1],
                    points_in_sphere[:,2],
                    c=activations2[distance<cutoff]-activations1[distance<cutoff],
                    marker='s',cmap='bwr',vmin=-3,vmax=3)
            fig.colorbar(p)
            ax.set_xlim3d(0,20)
            ax.set_ylim3d(0,20)
            ax.set_zlim3d(0,20)
            plt.title("Difference")
            plt.show()
            

    Relevant resources / publications

    Go back

    Multivoxel pattern classification

    description

    Upsides

    lorem

    Downsides

    Compared to IEMs, no tuning function no inherent difference betw classes Compared to NNs, given enough training data and computational power, NNs will always perform better due to universal approximation theorem

    Example

    First import the data and a few standard Python packages.
    
            # general packages
            import numpy as np
            import matplotlib.pyplot as plt
            import scipy as sp
            from sklearn import svm
            
            # pull data from GitHub
            import requests, io
            for array in ['activations','sub_id','cond_id','vox_id','block_id','trial_id']:
                globals()['{}'.format(array)] = np.load(io.BytesIO(requests.get(
                    'http://paulscotti.github.io/NavigateFMRI/methods/mvpa/mvpc/{}.npy'.format(array)).content))
            
            # overview of the data
            num_subjects = len(np.unique(sub_id)) #20
            num_conditions = len(np.unique(cond_id)) #2
            num_voxels = len(np.unique(vox_id)) #50
            num_blocks = len(np.unique(block_id)) #10
            
    We have 20 subjects who underwent 10 blocks, where each block was either watching a comedy movie clip or a horror movie clip. We have activations for each block for each of 50 voxels in our brain region of interest. Let's take a look at the overall histogram of activations for each condition (across all participants).
    
            plt.hist(activations[cond_id==0],color='blue',bins=30,label='Comedy')
            plt.hist(activations[cond_id==1],color='red',bins=30,alpha=.5,label='Horror')
            plt.yticks([]);plt.ylabel('Frequency')
            plt.xlabel("Activation strength")
            plt.legend(bbox_to_anchor=(1, 1))
            plt.show()
            

    It looks like there is a pretty clear difference between conditions, where the majority of comedy blocks were associated with higher activations. Let's train a model that can take activation and predict the presented condition (comedy or horror). To do this, we are going to use a linear support vector machine (SVM) to train on the even numbered blocks and test the accuracy of our model using the held-out, odd numbered blocks. SVMs are a kind of machine learning classifier where the objective of the classifier is to find a boundary (hyperplane) that maximizes the distance between classes (the margin). We will separately train and test the model for each subject, and then perform a t-test comparing all our prediction accuracies against what would be expected by chance.
    
            scores=[];
            for s in range(num_subjects):
                model = svm.SVC(kernel='linear')
            
                samples_train = activations[(sub_id==s) & (block_id%2==0)].reshape(-1, 1)
                classes_train = cond_id[(sub_id==s) & (block_id%2==0)].ravel()
            
                samples_test = activations[(sub_id==s) & (block_id%2==1)].reshape(-1, 1)
                classes_test = cond_id[(sub_id==s) & (block_id%2==1)].ravel()
            
                model.fit(samples_train,classes_train)
                scores = np.concatenate([scores,[model.score(samples_test,classes_test)]])
            plt.scatter(np.arange(len(scores)),scores)
            plt.ylim([.0,1.])
            plt.ylabel("% correct")
            plt.xlabel("Subject ID")
            plt.axhline(.5,ls='--',c='k',label='chance')
            plt.legend()
            plt.show()
            print("Avg correct predictions across subs:",np.mean(scores))
            stat = sp.stats.ttest_rel(scores,np.ones(num_subjects)*.5)
            print("Above-chance prediction accuracy? (t-test):\nt={:.4f}, p={:.4f}".format(stat.statistic,stat.pvalue))
            

    Avg correct predictions across subs: 0.8266
    Above-chance prediction accuracy? (t-test): t=9.0744, p=0.0000

    Let's visualize what one of our models is doing. For subject 1, we have 30 activations. 15 of these activations are used to train the model. The other 15 are held-out to test the model and are plotted as a histogram below. Our SVM uses the training data to estimate a hyperplane (aka decision boundary) that acts as a boundary between conditions (black vertical line). Activations to the right of the hyperplane are predicted to be one condition and activations to the left of the hyperplane are predicted to be the other condition. Note that this example is one-dimensional and is therefore easy to visualize, but hyperplanes can span dimensions and therefore SVMs can be used to classify much more complex data.
    
            s=0
            model = svm.SVC(kernel='linear')
            samples_train = activations[(sub_id==s) & (block_id%2==0)].reshape(-1, 1)
            classes_train = cond_id[(sub_id==s) & (block_id%2==0)].ravel()
            samples_test = activations[(sub_id==s) & (block_id%2==1)].reshape(-1, 1)
            classes_test = cond_id[(sub_id==s) & (block_id%2==1)].ravel()
            model.fit(samples_train,classes_train)
            
            plt.hist(activations[(sub_id==s) & (block_id%2==1) & (cond_id==0)],color='blue',bins=30,label='Comedy')
            plt.hist(activations[(sub_id==s) & (block_id%2==1) & (cond_id==1)],color='red',bins=30,alpha=.5,label='Horror')
            plt.yticks([]);plt.ylabel('Frequency')
            plt.xlabel("Activation strength")
            plt.legend(bbox_to_anchor=(1, 1))
            
            xx = np.linspace(-5, 5, 100)
            yy = np.linspace(0, 40, 100)
            YY, XX = np.meshgrid(yy, xx)
            xy = np.vstack([XX.ravel()]).T
            Z = model.decision_function(xy).reshape(XX.shape)
            
            plt.title("Subject 1")
            plt.contour(XX, YY, Z, colors='k', levels=[0],
                        linestyles=['--', '-', '--'])
            plt.show()
            

    Relevant resources / publications

  • Multi-Voxel Pattern Analysis of fMRI Data
  • Go back

    Multivoxel pattern analysis: correlation

    description

    Upsides

    lorem

    Downsides

    lorem

    Example

    First import the data and a few standard Python packages.
    
            # general packages
            import numpy as np
            import matplotlib.pyplot as plt
            import scipy as sp
            
            # pull data from GitHub
            import requests, io
            for array in ['activations','sub_id','cond_id','block_id','vox_id']:
                globals()['{}'.format(array)] = np.load(io.BytesIO(requests.get(
                    'http://paulscotti.github.io/NavigateFMRI/methods/mvpa/mvpcorrelation/{}.npy'.format(array)).content))
            
            # overview of the data
            num_subjects = len(np.unique(sub_id)) #20
            num_conditions = len(np.unique(cond_id)) #2
            num_voxels = len(np.unique(vox_id)) #50
            num_blocks = len(np.unique(block_id)) #40
            
    We have 20 subjects who underwent 40 blocks, where each block was either watching a comedy movie clip or a horror movie clip. We have activations for each block for each of 50 voxels in our brain region of interest. Let's visualize what this looks like for one subject.
    
            plt.imshow(activations[(sub_id==0)&(cond_id==0)].reshape(num_blocks//2,num_voxels),cmap='PuBu')
            plt.title("Comedy (Subject 1)")
            plt.ylabel("Block ID")
            plt.xlabel("Voxel ID")
            plt.colorbar()
            plt.show()
            plt.imshow(activations[(sub_id==0)&(cond_id==1)].reshape(num_blocks//2,num_voxels),cmap='OrRd')
            plt.title("Horror (Subject 1)")
            plt.ylabel("Block ID")
            plt.xlabel("Voxel ID")
            plt.colorbar()
            plt.show()
            

    The question we are asking with correlation-based MVPA is whether activations are more consistent within the same experimental condition compared to across different experimental conditions. First, we will divide the data in half, split by even or odd numbered blocks.
    
            plt.imshow(activations[(sub_id==0)&(cond_id==0)&(block_id%2==1)].reshape(num_blocks//4,num_voxels),cmap='gray')
            plt.title("Comedy (Sub 1, odd blocks)")
            plt.ylabel("Block")
            plt.xlabel("Voxel ID")
            plt.colorbar()
            plt.show()
            plt.imshow(activations[(sub_id==0)&(cond_id==0)&(block_id%2==0)].reshape(num_blocks//4,num_voxels),cmap='gray')
            plt.title("Comedy (Sub 1, even blocks)")
            plt.ylabel("Block")
            plt.xlabel("Voxel ID")
            plt.colorbar()
            plt.show()
            
            plt.imshow(activations[(sub_id==0)&(cond_id==1)&(block_id%2==1)].reshape(num_blocks//4,num_voxels),cmap='gray')
            plt.title("Horror (Sub 1, odd blocks)")
            plt.ylabel("Block")
            plt.xlabel("Voxel ID")
            plt.colorbar()
            plt.show()
            plt.imshow(activations[(sub_id==0)&(cond_id==1)&(block_id%2==0)].reshape(num_blocks//4,num_voxels),cmap='gray')
            plt.title("Horror (Sub 1, even blocks)")
            plt.ylabel("Block")
            plt.xlabel("Voxel ID")
            plt.colorbar()
            plt.show()
            


    Then, for each subject, we will measure the correlation in activations between even and odd blocks for every pairwise comparison of conditions (horror[odd] to horror[even], comedy[odd] to comedy[even], horror[odd] to comedy[even], and comedy[odd] to horror[even]).
    
            # for each subject, correlate odd/even runs for each condition
            corr = np.full((num_subjects,num_conditions*2),np.nan)
            for s in range(num_subjects):
                for c in range(num_conditions*2):
                    if c == 0: # Comedy (odd) -> Comedy (even)
                        corr[s,c], _ = sp.stats.pearsonr(activations[(sub_id==s)&(cond_id==0)&(block_id%2==1)],
                                                        activations[(sub_id==s)&(cond_id==0)&(block_id%2==0)])
                    elif c == 1: # Horror (odd) -> Horror (even)
                        corr[s,c], _ = sp.stats.pearsonr(activations[(sub_id==s)&(cond_id==1)&(block_id%2==1)],
                                                        activations[(sub_id==s)&(cond_id==1)&(block_id%2==0)])
                    elif c == 2: # Comedy (odd) -> Horror (even)
                        corr[s,c], _ = sp.stats.pearsonr(activations[(sub_id==s)&(cond_id==0)&(block_id%2==1)],
                                                        activations[(sub_id==s)&(cond_id==1)&(block_id%2==0)])
                    elif c == 3: # Horror (odd) -> Comedy (even)
                        corr[s,c], _ = sp.stats.pearsonr(activations[(sub_id==s)&(cond_id==1)&(block_id%2==1)],
                                                        activations[(sub_id==s)&(cond_id==0)&(block_id%2==0)])
                
            with plt.rc_context(rc={'font.size': 20, 'figure.figsize': (5,2)}):
                plt.imshow(corr, cmap='gray', aspect='auto')
                plt.ylabel('Subject ID')
                plt.yticks([]);
                plt.xticks(np.arange(4),["Horror\n(Odd)\nHorror\n(Even)","Comedy\n(Odd)\nComedy\n(Even)",
                                            "Horror\n(Odd)\nComedy\n(Even)","Comedy\n(Odd)\nHorror\n(Even)"],
                            fontsize=14, rotation=0, ha="center")
                plt.colorbar()
                plt.title("Correlation coefficients")
                plt.show()
            
            print('From left to right:')
            print('° Comedy (odd) & Comedy (even)  Avg. across subjects:  r={:.4f}'.format(np.mean(corr[:,0])))
            print('° Horror (odd) & Horror (even)  Avg. across subjects:  r={:.4f}'.format(np.mean(corr[:,1])))
            print('° Comedy (odd) & Horror (even)  Avg. across subjects:  r={:.4f}'.format(np.mean(corr[:,2])))
            print('° Horror (odd) & Comedy (even)  Avg. across subjects:  r={:.4f}'.format(np.mean(corr[:,3])))
            
            stat = sp.stats.ttest_rel(corr[:,:2].flatten(),corr[:,2:].flatten())
            print("\nSignificantly different correlations betw first two columns and last two columns? (paired t-test):\nt={:.4f}, p={:.4f}".format(stat.statistic,stat.pvalue))
            

  • Comedy (odd) & Comedy (even): Avg. across subjects: r=0.2428
  • Horror (odd) & Horror (even): Avg. across subjects: r=0.3411
  • Comedy (odd) & Horror (even): Avg. across subjects: r=-0.2605
  • Horror (odd) & Comedy (even): Avg. across subjects: r=-0.2612
  • Significantly different correlations betw first two columns and last two columns? (paired t-test): t=48.5594, p=0.0000

    As we can see in the above figure, the correlations between same-condition blocks are much higher than the correlations between different-condition blocks. A t-test shows that these differences are significant. That is, the pattern of activations in this brain region are showing differences dependent on the experimental condition. Therefore, there is some condition-specific processing occurring in this brain region (this brain region responds differently to comedy movies than horror movies).

    Note that in this example we have not actually classified any data (as we do in multivariate pattern classification), but it is possible to use correlation-based MVPA for classification. For instance, you could split your data in half, where in one half you know the labels for your data (training set) and in the other half you do not know the labels of your data (test set). And then, you can correlate your test set with your training set and predict the labels of the test set based on the highest correlation with the training set.

    Also note pairwise correlations of conditions are often depicted in what is known as a representational similarity matrix (RSM). This becomes especially relevant for Representational Similarity Analysis, which involves comparing RSMs. Below is the RSM for subject 1 between odd and even numbered blocks.
    
                sub=0 # for first subject
                # Calcuate correlations between even- and odd-block conditions
                comedy_even = activations[(sub_id==0)&(cond_id==0)&(block_id%2==0)]
                comedy_odd = activations[(sub_id==0)&(cond_id==0)&(block_id%2==1)]
                horror_even = activations[(sub_id==0)&(cond_id==1)&(block_id%2==0)]
                horror_odd = activations[(sub_id==0)&(cond_id==1)&(block_id%2==1)]
                # Plot correlations between even- and odd-block conditions
                correlation_table = []
                # first row of correlation matrix:
                correlation_table.append([stats.pearsonr(comedy_even,comedy_odd)[0],
                                        stats.pearsonr(horror_even,comedy_odd)[0]])
                # second row of correlation matrix:
                correlation_table.append([stats.pearsonr(comedy_even,horror_odd)[0],
                                        stats.pearsonr(horror_even,horror_odd)[0]])
    
                # plot correlation matrix
                plt.rcParams.update({'font.size': 20, 'figure.figsize': (5,5)}) # change default plotting
                plt.imshow(correlation_table,cmap='hot')
                plt.title("Representational Similarity Matrix\n(Subject {})".format(sub+1),fontsize=18)
                plt.yticks(np.arange(2),["Comedy","Horror"],fontsize=16)
                plt.ylabel("Odd")
                plt.xticks(np.arange(2),["Comedy","Horror"],fontsize=16)
                plt.xlabel("Even")
                cbar = plt.colorbar()
                cbar.set_label('Pearson R', rotation=270)
                plt.show()
            

    Relevant resources / publications

  • Multi-Voxel Pattern Analysis of fMRI Data
  • Go back

    Representational similarity analysis

    description

    Upsides

    lorem

    Downsides

    lorem

    Example

    First import the data and a few standard Python packages.
    
                # general packages
                import numpy as np
                import matplotlib.pyplot as plt
                import scipy as sp
    
                # pull data from GitHub
                import requests, io
                for array in ['activations','sub_id','cond','block_id','actor']:
                    globals()['{}'.format(array)] = np.load(io.BytesIO(requests.get(
                        'http://paulscotti.github.io/navigateFMRI/methods/mvpa/rsa/{}.npy'.format(array)).content))
            
    We have 20 subjects each with a 50-voxel brain region of interest (ROI). We have activations for every voxel for each of four conditions (Comedy starring Jennifer Aniston, Comedy starring Johnny Depp, Horror starring Jennifer Aniston, and Horror starring Johnny Depp).

    As briefly introduced in the correlation-based MVPA example, RSA involves creating a representational similarity matrix (RSM). An RSM is simply a matrix of all pairwise correlations (or any sort of distance measure) of conditions (e.g., splitting your data in half and correlating each half per condition).
    
                sub=0 # for first subject
                # Calcuate correlations between even- and odd-block conditions
                comedy_aniston_even = activations[(sub_id==sub)&(cond=='comedy')&(block_id%2==0)&(actor=="Aniston")]
                comedy_aniston_odd = activations[(sub_id==sub)&(cond=='comedy')&(block_id%2==1)&(actor=="Aniston")]
                comedy_depp_even = activations[(sub_id==sub)&(cond=='comedy')&(block_id%2==0)&(actor=="Depp")]
                comedy_depp_odd = activations[(sub_id==sub)&(cond=='comedy')&(block_id%2==1)&(actor=="Depp")]
                horror_aniston_even = activations[(sub_id==sub)&(cond=='horror')&(block_id%2==0)&(actor=="Aniston")]
                horror_aniston_odd = activations[(sub_id==sub)&(cond=='horror')&(block_id%2==1)&(actor=="Aniston")]
                horror_depp_even = activations[(sub_id==sub)&(cond=='horror')&(block_id%2==0)&(actor=="Depp")]
                horror_depp_odd = activations[(sub_id==sub)&(cond=='horror')&(block_id%2==1)&(actor=="Depp")]
    
                # Plot correlations between even- and odd-block conditions
                correlation_table = []
                # first row of correlation matrix:
                correlation_table.append([stats.pearsonr(comedy_aniston_even,comedy_aniston_odd)[0],
                                          stats.pearsonr(comedy_depp_even,comedy_aniston_odd)[0],
                                          stats.pearsonr(horror_aniston_even,comedy_aniston_odd)[0],
                                          stats.pearsonr(horror_depp_even,comedy_aniston_odd)[0]])
                # second row of correlation matrix:
                correlation_table.append([stats.pearsonr(comedy_aniston_even,comedy_depp_odd)[0],
                                          stats.pearsonr(comedy_depp_even,comedy_depp_odd)[0],
                                          stats.pearsonr(horror_aniston_even,comedy_depp_odd)[0],
                                          stats.pearsonr(horror_depp_even,comedy_depp_odd)[0]])
                # third row of correlation matrix:
                correlation_table.append([stats.pearsonr(comedy_aniston_even,horror_aniston_odd)[0],
                                          stats.pearsonr(comedy_depp_even,horror_aniston_odd)[0],
                                          stats.pearsonr(horror_aniston_even,horror_aniston_odd)[0],
                                          stats.pearsonr(horror_depp_even,horror_aniston_odd)[0]])
                # fourth row of correlation matrix:
                correlation_table.append([stats.pearsonr(comedy_aniston_even,horror_depp_odd)[0],
                                          stats.pearsonr(comedy_depp_even,horror_depp_odd)[0],
                                          stats.pearsonr(horror_aniston_even,horror_depp_odd)[0],
                                          stats.pearsonr(horror_depp_even,horror_depp_odd)[0]])
                # plot correlation matrix
                plt.rcParams.update({'font.size': 20, 'figure.figsize': (5,5)}) # change default plotting
                plt.imshow(correlation_table,cmap='hot')
                plt.title("Representational Similarity Matrix\n(Subject {})".format(sub+1))
                plt.yticks(np.arange(4),["Aniston","Depp","Aniston","Depp"],fontsize=16)
                plt.ylabel("Horror-Odd  Comedy-Odd")
                plt.xticks(np.arange(4),["Aniston","Depp","Aniston","Depp"],fontsize=16)
                plt.xlabel("Comedy-Even  Horror-Even")
                cbar = plt.colorbar()
                cbar.set_label('Pearson R', rotation=270)
                plt.show()
            

    Let us say that we have two theoretical questions: (1) The ROI represents movie genre; (2) The ROI represents specific actors. These theoretical questions can be represented as RSMs like so:
    
                ComedyVsHorror= []
                ComedyVsHorror.append([1,1,0,0])
                ComedyVsHorror.append([1,1,0,0])
                ComedyVsHorror.append([0,0,1,1])
                ComedyVsHorror.append([0,0,1,1])
                # plot correlation matrix
                plt.rcParams.update({'font.size': 20, 'figure.figsize': (5,5)}) # change default plotting
                plt.imshow(ComedyVsHorror,cmap='RdYlBu')
                plt.title("Genre RSM")
                plt.yticks(np.arange(4),["Aniston","Depp","Aniston","Depp"],fontsize=16)
                plt.ylabel("Horror        Comedy")
                plt.xticks(np.arange(4),["Aniston","Depp","Aniston","Depp"],fontsize=16)
                plt.xlabel("Comedy        Horror")
                plt.show()
                
                AnistonVsDepp = []
                AnistonVsDepp.append([1,0,1,0])
                AnistonVsDepp.append([0,1,0,1])
                AnistonVsDepp.append([1,0,1,0])
                AnistonVsDepp.append([0,1,0,1])
                # plot correlation matrix
                plt.rcParams.update({'font.size': 20, 'figure.figsize': (5,5)}) # change default plotting
                plt.imshow(AnistonVsDepp,cmap='RdYlBu')
                plt.title("Actor RSM")
                plt.yticks(np.arange(4),["Aniston","Depp","Aniston","Depp"],fontsize=16)
                plt.ylabel("Horror        Comedy")
                plt.xticks(np.arange(4),["Aniston","Depp","Aniston","Depp"],fontsize=16)
                plt.xlabel("Comedy        Horror")
                plt.show()
            

    We can then test how similar these RSMs are to our obtained, simulated data. Essentially, this is just performing a t-test between the blue and red cells for each of the two research questions. This is called 2nd-order RSA (aka representational connectivity, or informational connectivity).
    
                # Calcuate correlations between even- and odd-block conditions
                comedy_allsub = np.full((num_subjects),np.nan)
                horror_allsub = np.full((num_subjects),np.nan)
                aniston_allsub = np.full((num_subjects),np.nan)
                depp_allsub = np.full((num_subjects),np.nan)
                for sub in range(num_subjects):
                    correlation_table=[]
                    comedy_aniston_even = activations[(sub_id==sub)&(cond=='comedy')&(block_id%2==0)&(actor=="Aniston")]
                    comedy_aniston_odd = activations[(sub_id==sub)&(cond=='comedy')&(block_id%2==1)&(actor=="Aniston")]
                    comedy_depp_even = activations[(sub_id==sub)&(cond=='comedy')&(block_id%2==0)&(actor=="Depp")]
                    comedy_depp_odd = activations[(sub_id==sub)&(cond=='comedy')&(block_id%2==1)&(actor=="Depp")]
                    horror_aniston_even = activations[(sub_id==sub)&(cond=='horror')&(block_id%2==0)&(actor=="Aniston")]
                    horror_aniston_odd = activations[(sub_id==sub)&(cond=='horror')&(block_id%2==1)&(actor=="Aniston")]
                    horror_depp_even = activations[(sub_id==sub)&(cond=='horror')&(block_id%2==0)&(actor=="Depp")]
                    horror_depp_odd = activations[(sub_id==sub)&(cond=='horror')&(block_id%2==1)&(actor=="Depp")]
                    # first row of correlation matrix:
                    correlation_table.append([sp.stats.pearsonr(comedy_aniston_even,comedy_aniston_odd)[0],
                                              sp.stats.pearsonr(comedy_depp_even,comedy_aniston_odd)[0],
                                              sp.stats.pearsonr(horror_aniston_even,comedy_aniston_odd)[0],
                                              sp.stats.pearsonr(horror_depp_even,comedy_aniston_odd)[0]])
                    # second row of correlation matrix:
                    correlation_table.append([sp.stats.pearsonr(comedy_aniston_even,comedy_depp_odd)[0],
                                              sp.stats.pearsonr(comedy_depp_even,comedy_depp_odd)[0],
                                              sp.stats.pearsonr(horror_aniston_even,comedy_depp_odd)[0],
                                              sp.stats.pearsonr(horror_depp_even,comedy_depp_odd)[0]])
                    # third row of correlation matrix:
                    correlation_table.append([sp.stats.pearsonr(comedy_aniston_even,horror_aniston_odd)[0],
                                              sp.stats.pearsonr(comedy_depp_even,horror_aniston_odd)[0],
                                              sp.stats.pearsonr(horror_aniston_even,horror_aniston_odd)[0],
                                              sp.stats.pearsonr(horror_depp_even,horror_aniston_odd)[0]])
                    # fourth row of correlation matrix:
                    correlation_table.append([sp.stats.pearsonr(comedy_aniston_even,horror_depp_odd)[0],
                                              sp.stats.pearsonr(comedy_depp_even,horror_depp_odd)[0],
                                              sp.stats.pearsonr(horror_aniston_even,horror_depp_odd)[0],
                                              sp.stats.pearsonr(horror_depp_even,horror_depp_odd)[0]])
                    comedy_allsub[sub] = np.mean(np.array(correlation_table)[np.array(ComedyVsHorror)==1])
                    horror_allsub[sub] = np.mean(np.array(correlation_table)[np.array(ComedyVsHorror)==0])
                    aniston_allsub[sub] = np.mean(np.array(correlation_table)[np.array(AnistonVsDepp)==1])
                    depp_allsub[sub] = np.mean(np.array(correlation_table)[np.array(AnistonVsDepp)==0])
                    
                plt.bar(0,np.mean(comedy_allsub),color='blue')
                plt.bar(1,np.mean(horror_allsub),color='darkred')
                plt.scatter(np.zeros(num_subjects),comedy_allsub,color='k',zorder=2)
                plt.scatter(np.ones(num_subjects),horror_allsub,color='k',zorder=2)
                plt.axhline(0,c='k')
                plt.title("Same-genre Vs Different-genre")
                plt.xticks(np.arange(2),["Same","Diff."])
                plt.ylabel("Avg. correlation")
                plt.ylim(0,.1)
                plt.show()
                stat = sp.stats.ttest_rel(comedy_allsub,horror_allsub)
                print("\nPaired t-test: t={:.4f}, p={:.4f}".format(stat.statistic,stat.pvalue))
                
                plt.bar(0,np.mean(aniston_allsub),color='blue')
                plt.bar(1,np.mean(depp_allsub),color='darkred')
                plt.scatter(np.zeros(num_subjects),aniston_allsub,color='k',zorder=2)
                plt.scatter(np.ones(num_subjects),depp_allsub,color='k',zorder=2)
                plt.axhline(0,c='k')
                plt.title("Same-actor Vs Different-actor")
                plt.xticks(np.arange(2),["Same","Diff."])
                plt.ylabel("Avg. correlation")
                plt.ylim(0,.1)
                plt.show()
                stat = sp.stats.ttest_rel(aniston_allsub,depp_allsub)
                print("\nPaired t-test: t={:.4f}, p={:.4f}".format(stat.statistic,stat.pvalue))
            

    Paired t-test: t=5.5199, p=0.0000 Paired t-test: t=-0.0462, p=0.9636

    Based on these results, we can conclude that the ROI did differentially represent movie genre but there was not sufficient evidence to suggest differentially represented actors.

    Relevant resources / publications

  • Dr. Rebecca Saxe fMRI Bootcamp Part 7
  • Kriegeskorte, Mur, & Bandettini, 2008
  • Go back

    Functional connectivity

    description

    Upsides

    lorem

    Downsides

    lorem

    Example

    First import the data and a few standard Python packages.
    
                # general packages
                import numpy as np
                import matplotlib.pyplot as plt
                import scipy as sp
                
                # pull data from GitHub
                import requests, io
                for array in ['enjoyment','time_series1','time_series2']:
                    globals()['{}'.format(array)] = np.load(io.BytesIO(requests.get(
                        'http://paulscotti.github.io/navigateFMRI/methods/univariate/{}.npy'.format(array)).content))
                
                # overview of the data
                num_subjects = time_series1.shape[0]
                num_TRs = time_series1.shape[1]
                num_conditions = 2
            
    In this simulated dataset we have 20 subjects. Each subject has a BOLD signal time series for two brain regions (ROI A and ROI B) which represents the brain activations for each subject watching a movie. Each subject also rated their enjoyment of the movie on a scale from 1 to 5. One question we can ask is whether stronger functional connectivity between the two brain regions varies with participants' subjective enjoyment of the movie. First let's plot the time series of each brain region for one subject, where the x-axis here represents the repetition time (TR) of the scanner (MRI scanners only collect data every x seconds, where x is the repetition time).
    
                plt.plot(time_series1[0,:],c='green')
                plt.title("Subject 1, ROI A")
                plt.ylabel("BOLD signal")
                plt.xlabel("TR")
                plt.show()
                plt.plot(time_series2[0,:],c='purple')
                plt.title("Subject 1, ROI B")
                plt.ylabel("BOLD signal")
                plt.xlabel("TR")
                plt.show()
                print("Correlation coefficient between the two time series = {}".format(
                        sp.stats.pearsonr(time_series1[0,:],time_series2[0,:])[0]))
                print("Movie enjoyment rating = {} (1-5 scale)".format(enjoyment[0]))
            

    Correlation coefficient between the two time series = 0.1789
    Movie enjoyment rating = 4.0 (1-5) scale

    We can see that there is a positive correlation between the two time-series, but is this relevant for behavior? Let's plot the correlation coefficient for every subject along with the movie enjoyment rating for every subject.
    
                corrs=[]
                for sub in range(num_subjects):
                    coefficient=sp.stats.pearsonr(time_series1[sub,:],time_series2[sub,:])[0]
                    plt.scatter(sub,coefficient,color='brown')
                    corrs=np.append(corrs,coefficient)
                plt.title("Time-series correlation coefficient")
                plt.xlabel("Subject")
                plt.show()
                for sub in range(num_subjects):
                    plt.scatter(sub,enjoyment[sub],color='black')
                plt.title("Movie enjoyment rating")
                plt.yticks(np.arange(5)+1,fontsize=14)
                plt.xlabel("Subject")
                plt.show()
            

    Visually it seems there might be some association here, where higher correlation coefficients are associated with increased movie enjoyment. Finally, let's test if this is statistically significant by using a paired t-test.
    
                # do statistics between conditions
                stat = sp.stats.ttest_rel(corrs,enjoyment)
                print("Correlation coefficients sig. associated with enjoyment ratings? (paired t-test):\nt={:.4f}, p={:.4f}".format(stat.statistic,stat.pvalue))
            
    Correlation coefficients sig. associated with enjoyment ratings? (paired t-test):
    t=-11.0475, p=0.0000

    So it seems that increased functional connectivity is correlated with movie enjoyment. But be careful when interpreting such results -- just because a correlation is significant does not imply anything about the cause or directionality of the effect. For instance, you cannot conclude that increased connectivity between these two brain regions causes increased movie enjoyment, or even that these brain regions are involved at all in the processing of movie enjoyment. It's possible that movie enjoyment is correlated with increased attention, for instance, and that attention leads to downstream increased processing within these brain regions.

    Relevant resources / publications

  • Rogers, B. P., Morgan, V. L., Newton, A. T., & Gore, J. C. (2007). Assessing functional connectivity in the human brain by fMRI. Magnetic resonance imaging, 25(10), 1347-1357.
  • Hampson M, Tokoglu F, Sun Z, Schafer RJ, Skudlarski P, Gore JC, Constable RT. Connectivity-behavior analysis reveals that functional connectivity between left BA39 and Broca's area varies with reading ability. Neuroimage. 2006; 31(2), 513-519.
  • Go back

    Method name

    description

    Upsides

    lorem

    Downsides

    lorem

    Example

    First import the data and a few standard Python packages.
    
                code
            

    Relevant resources / publications

    Go back

  • Regarding general linear models, general refers to the different types of analyses that can be used (e.g., correlation, one-sample t-test, ANOVA, etc.) and linear refers to how the model assumes a linear relationship between parameters (often preferred because nonlinear models are prone to overfitting, where a model appears to fit the data well but is not generalizable to new data). We consider the general linear model to be a processing step (rather than a post-processing fMRI method, which is what this website is focused on) and thus will not discuss it in further detail, but it is essentially one way to obtain an estimate of brain activation (in arbitrary units) for every voxel for each condition after accounting for things like the hemodynamic response function (time-course of blood flow) and body motion in the scanner.

  • In practice, you should use cross-validation such that all data may be predicted without succumbing to circularity/double-dipping.

  • Unlike NavigateFMRI, Brainiak involves installing a custom Python package and is built to be a helpful tool for use on real data. NavigateFMRI is meant purely for teaching purposes with no intention for the code contained on the website to be directly implemented in real fMRI projects. Brainiak's Python tutorial examples use real fMRI datasets and go into more depth than NavigateFMRI. NavigateFMRI covers more methods than Brainiak and uses simpler, simulated datasets. Our example walkthroughs are much shorter and do not dwell on technical details that may be important to consider when actually implementing techniques on real data.

  • NavigateFMRI only covers post-processing fMRI methods while fmri4newbies covers all stages of the fMRI pipeline including pre-processing, processing, and post-processing (in addition to information on MRI physics and anatomical & functional brain regions). Unlike NavigateFMRI, fmri4newbies supplies the reader with Powerpoint lecture slides (no example datasets or code).

  • "Activations" can be raw BOLD signal or beta weights obtained from general linear model estimation.

  • Please just pretend that Jennifer Aniston has starred in horror movies, we are not movie buffs :). Note: Jennifer Aniston actually has a reputation for being used as an example in neuroimaging experiments because of a study investigating "Grandmother cells", the idea that a neuron may be so complex as to specifically encode your grandmother. The study involved one participant who seemed have a neuron with Jennifer Aniston specific firing.