RSA: how to describe with a single number? - update

Some years ago I posted about quantifying RSA matrices with a single number. Since then, it seems that the field has settled on two methods for quantifying RSA matrices: one based on differences and the other based on Kendall's tau-a.

To illustrate the methods, suppose this is the RSA matrix for one person, and the corresponding "reference matrix"; the numbers in the cells are Pearson correlations. In this case, we expect the two cells with the same letter (2,F and 0,F; 2,P and 0,P) to have higher correlations than all the other cells (which mix letters).

The matrix is symmetrical, with 1s along the diagonal, so we'll only use the numbers from the lower triangle. We can "unwrap" these values into vectors, like this:

 x <- c(.728, .777, .705, .709, .873, .705);   # one person
 y <- c(   0,    1,    0,    0,    1,    0);   # reference

I previously called the difference-based method "averaging the triangles"; conceptually, we take the difference between the average of the two cells we expect to have a large value (cells with 1s in the Reference), and the average of the four cells we expect to have smaller values (0 in the Reference). This can also be thought of as using a contrast matrix, such as in Oedekoven, et. al (2017). Since you can't do math on correlations without going through the Fisher R-to-Z transform, here's an example of the calculation:

 get.FTrz <- function(in.val) { return(.5 * log((1+in.val)/(1-in.val))); }  # Fisher's r-to-z transformation  
 get.FTzr <- function(in.val) { return((exp(2*in.val)-1)/(exp(2*in.val)+1)); } # Fisher's z-to-r transformation  
 mean(get.FTrz(x[which(y == 1)])) - mean(get.FTrz(x[which(y == 0)]));
 # [1] 0.3006614  

(It's of course equivalent to calculate this by multiplying the sum of the y == 1 cells by 0.5 and y == 0 cells by 0.25.)

The Kendall's tau-a based method is recommended in Nili, et. al (2014). There's a function to calculate this in the DescTools R package; since it's rank-based, there's no need to do the Fisher transformation:

 KendallTauA(x,y)  
 # [1] 0.5333333  

Which method should you use? One consideration is if you want the magnitude of the numbers in the RSA matrix cells to matter, or just their relative sizes (ranks). Here, I made a RSA matrix with the same rank ordering as before (in x), but the difference between same-letter and different-letter cells are more extreme:

 x1 <- c(0.2, 0.9, 0.1, 0.01, 0.83, 0.22);  
 KendallTauA(x1,y)  
 # [1] 0.5333333  # same  
 mean(get.FTrz(x1[which(y == 1)])) - mean(get.FTrz(x1[which(y == 0)]));
 # [1] 1.195997 # bigger  

The Kendall's tau-a value for x1 and y is exactly the same as x and y, but the mean difference is larger. Is it better to be sensitive to magnitude? Not necessarily. In practice, when I've calculated both methods on the same dataset I've found extremely similar results, which seems proper; it's hard to imagine a strong effect that would be only be detected by Kendall's tau-a.

UPDATE 12 December 2018: I removed transforming the difference scores back to correlations, leaving the differences as Fisher-transformed zs, since this seems easier for later interpretation and calculations.

detrending and normalizing timecourses: afni 3dDetrend in R

I've been using the afni 3dDetrend command on datasets lately, but realized that I wasn't sure of the exact way the normalization and detrending was calculated. To help me understand (and reproduce the results on files that aren't easily read by afni), this post gives R code to match 3dDetrend output, plus shows some examples of how various normalizing and detrending schemes change timecourses (there's more in the knitr and pdf than shown in this post; files below the jump).

First, we need some data. I picked a run from a handy preprocessed dataset (it doesn't really matter for the demo, but it's from multibandACQtests), then ran the image through afni's 3dDetrend command twice: once with the options -normalize -polort 0, then with -normalize -polort 2. This gives three images, and I picked a voxel at random and extracted its timecourse from each image. Code for this step starts at line 22 of the knitr file (below the jump), and text files with the values are included in the archive for this post.

Here is the voxel's "raw" timecourse: after preprocessing (the HCP pipelines, in this case), but before running 3dDetrend:

And here is the same voxel's timecourse, from the image after 3dDetrend -normalize -polort 2: centered on zero, with the shift of baseline around TR 50 reduced a bit.

From the afni documentation, 3dDetrend -normalize "Normalize[s] each output voxel time series; that is, make the sum-of-squares equal to 1." Note that this is not what R scale() does by default (which is to mean 0 and standard deviation 1; see the full demo for a comparison). Line 113 of the demo code does the afni-style normalization: (tc.raw - mean(tc.raw)) / sqrt(sum((tc.raw - mean(tc.raw))^2)), where tc.raw is a vector with the voxel's timecourse.

For the -polort 2 part, the afni documentation says: "Add Legendre polynomials of order up to and including 'ppp' in the list of vectors to remove." I used the functions from Gregor Gorjanc's blog post (thanks!) and the R orthopolynom package for this; the key part is line 140: residuals(lm(tc.raw ~ seq(tc.raw) + Legendre(x=seq(tc.raw), n=1) + Legendre(x=seq(tc.raw), n=2))). The two Legendre function terms (n=1 and n=2) are to match -polort 2; more (or fewer) terms can be included as desired.

tutorial: plotting GIfTI images in R (knitr)

NOTICE: this post was replaced by a newer version 13 March 2020. I will leave this post here, but strongly suggest you follow the new version instead.

In a previous tutorial I showed how to plot volumetric NIfTI images in R and knitr. Here, I give functions for plotting surface GIfTI brain images in R. If you're not familiar with surface images, this and this post may be useful. Please cite if you find this useful.

This uses the R gifti package (thank you, John Muschelli!) to read the GIfTI images, and the plot3D package for the plotting. It is possible to plot interactive 3d images in R (see for example, the gifti package vignette), but static images are more useful for the application I have in mind: batch-produced knitr files of surface (afni) GLM results.

Here is the first part of the knitr file produced in this tutorial; the full pdf can be downloaded here:

The "blobs" are F-statistics from a GLM targeting right-hand button pushes, so the peaks in left motor areas are sensible. The underlay surfaces are the "inflated" images from the HCP 1200 Subjects Group Average Data release; download here. The code requires you to have underlay surfaces and statistics to overlay in perfect agreement (the same number of nodes, same triangle connections, etc.), but this is a general requirement for surface plotting or analysis.

The gifti-plotting function I wrote displays the medial and lateral views, with hot colors for positive values and cool colors for negative. The function takes the positive and/or negative color scaling limits, with values closer to zero than the minimum value not shown, but those above the maximum plotted in the maximum color. This behavior is conventional for neuroimaging; for example, this is the same image, plotted in the Connectome Workbench with the same color scaling.

Full code for the demo is below the jump (or downloaded here). To run the demo you will need the two gifti overlay images, which can be downloaded here (left) and here (right), as well as the HCP 1200 underlay anatomies (see above). If you don't want to run the code within knitr the relevant functions can be separated from the latex sections, but you will need to set the image sizing and resolution.

Please let me know if you find this code useful, particularly if you develop any improvements!

UPDATE 20 July 2018: Fixed the plot.surface function to display the entire range (not truncated) in the lower right note.

UPDATE 7 September 2018: Another fix to plot.surface() to truncate the heat.colors() range so that the top values are plotted in yellow, rather than white, which wasn't very visible.

Also, Chris_Rorden posted a python and Surfice version on Neurostars.

UPDATE 9 November 2018: Added a section assigning arbitrary values to surface parcels and then plotting, matching this demo.

UPDATE 9 January 2019: Added a which.surface term to the plot.surface function, specifying whether plotting HCP or fsaverage5 surfaces. The default is HCP, for backward compatibility.

UPDATE 10 April 2019: Improved color ranges.

NOTICE: this post was replaced by a newer version 13 March 2020. I will leave this post here, but strongly suggest you follow the new version instead.

mean pattern subtraction confusion

A grad student brought several papers warning against subtracting the mean pattern (and other types of normalization) before correlational analyses (including, but not only, RSA) to my attention. I was puzzled: Pearson correlation ignores additive and  multiplicative transformations, so how could it be affected by subtracting the mean? Reading closely, my confusion was from what exactly was meant by "subtracting the mean pattern"; it is still the case that not all forms of "normalization" before correlational MVPA are problematic.


The key is where and how the mean-subtraction and/or normalization is done. Using the row- and column-scaling terminology from Pereira, Mitchell, and Botvinick (2009) (datasets arranged with voxels in columns, examples (trials, frames, whatever) in rows): the mean pattern subtraction warned against in Walther et al. (2016) and Garrido et al. (2013) is a particular form of column-scaling, not row-scaling. (A different presentation of some of these same ideas is in Hebart and Baker (2018), Figure 5.)

Here's my illustration of two different types of subtracting the mean; code to generate these figures is below the jump. The original pane shows the "activity patterns" for each of five trials in a four-voxel ROI. The appearance of these five trials after each has been individual normalized (row-scaled) is in the "trial normalized" pane, while the appearance after the mean pattern was subtracted is at right (c.f. Figure 1D in Walther et al. (2016) and Figure 1 in Garrido et al. (2013)).

Note that Trial2 is Trial1+15; Trial3 is Trial1*2; Trial4 is the mirror image of Trial1; Trial5 has a shape similar to Trial4 but positioned near Trial1. (Example adapted from Figure 8.1 of Cluster analysis for researchers (1984) by H. Charles Romesburg.)

And here are matrix views of the Pearson correlation and Euclidean distance between each pair of trials in each version of the dataset:

The mean pattern subtraction did not change the patterns' appearance as much as the trial normalization did: the patterns are centered on 0, but the extents are the same (e.g., voxel2 ranges from 0 to 80 in the original dataset, -40 to 40 after mean pattern subtraction). The row (trial) normalized image has a much different appearance: centered on zero, but the patterns for trials 1, 2, and 3 are now identical, and the range is much less (e.g., voxel2 is now a bit less than -1 to a bit more than 1). Accordingly, the trial-normalized similarity matrix is identical to the original when calculated with correlation, but different when calculated with Euclidean distance; the mean pattern subtracted similarity matrix is identical to the original when calculated with Euclidean distance, but different when calculated with Pearson correlation.

This is another case where lack of clear terminology can lead to confusion; "mean pattern subtraction" and "subtracting the mean from each pattern" are quite different procedures and have quite different effects, depending on your distance metric (correlation or Euclidean).

more crescents: with sequence variations

Back in January I posted about the "crescent" artifacts that show up in the PA runs of some people. (In a few people reversed crescents appear in the AP runs as well.) The working hypothesis is still that these are N/2 ghosts and perhaps related to insufficient fat suppression. I am still extremely interested in any advice people have for avoiding these, especially as I now have evidence that task signal is noticeably reduced in the "crescent" areas (details soon, hopefully).

We have begun a set of pilot scans to see if some parameter combinations produce better subcortical and frontal signal in a reward task. So far all scans have been on a Siemens Prisma, 64 32 channel head coil, CMRR MB4 sequences; we'll be doing the tests on a Siemens Vida as well in a few weeks. So far, scans are acquired with 2.4 or 3.0 mm isotropic voxels, either "flat" (AC-PC aligned as usual) or "tilted" (30 degrees off AC-PC); more acquisition details below the jump.

Of the two pilot people (so far), one has the crescent artifact and the other does not, and the appearance of the crescents in the different acquisitions is interesting. All of these images are voxelwise standard deviation, calculated over the entire run (no censoring, but extremely low motion), and on raw images (preprocessing is in progress).

First, here are sagittal views of a flat (left, scan 15) and tilted (right, scan 37) 2.4 mm isotropic run. The crescent artifact is visible in both; I marked the approximate ends with green arrows (click to enlarge). The multiband slice boundary is visible in both a fourth of the way up the image (red arrows).

Here are axial slices of the set of runs we have so far for this person. All are with the same color scaling (0 to 200); brighter is higher standard deviation. These are raw images, so the slice appearance varies quite a bit between the "flat" and "tilted" runs.

The "crescents" are visible in all PA runs, though perhaps easiest to spot with the 2p4 (2.4 mm isotropic) voxels. The slices in which the crescents appear varies between the tilted and flat acquisitions (e.g., k=31-46 for run 15_2p4flat_PA; 22-31 for run 37_2p4tilt_PA). It will be easier to compare the crescent locations after preprocessing.

The 3p0 (3.0 mm isotropic voxels) images are generally more uniform and dark than the 2p4 runs, likely reflecting improved signal-to-noise, particularly in the middle of the brain. While the large vessels are brightest in all runs (as they should be), the runs with 2.4 mm voxels (2p4) generally have a "starburst" type effect (brighter in the center, darker towards the edges), which is worrying, particularly since we want good signal in reward areas.

I will share other observations on this blog as the piloting and analyses progress. Please contact me if you'd like to run your own analyses; we'd be happy to share and are very interested in others' thoughts.

UPDATE 18 April 2018: I've wondered before if head size was a factor in which people have the crescent artifact, using the total intracranial volume measurement produced by freesurfer as a proxy. I don't have those measurements yet, but they kindly allowed me to measure their heads as if fitting them for hats, and they were nearly identical: about 58 cm for the person without the crescent artifact, and about 57 for the person shown in this person (with the artifact). Both people have a normal healthy body size; the person without the artifact was a bit shorter (around 5'2") than the other (around 5'7"). So, at least for these two pilots, external head size doesn't seem to matter for the artifact.

more acquisition parameters below the jump

Connectome Workbench: making a surface version of a volumetric image

This is a tutorial for making a surface version of a volumetric (NIfTI) image (e.g., a ROI) for visualizing with the Connectome Workbench. This post replaces a tutorial I wrote back in 2013, and assumes a bit of familiarity with Workbench. If you're using Workbench for the first time I suggest that you complete the tutorial in this post before trying this one.

First, a warning: as advertised, the steps in this post will make a surface version of a volumetric image. However, the surface version will be an approximation, and likely only suitable for visualization purposes (e.g., making an illustration of a set of ROIs for a talk). If you have an entire dataset that you want to prepare for surface analysis (e.g., running GLMs), you need different procedures (e.g., SUMA, FreeSurfer). Again, I suggest the directions in this post (wb_command -volume-to-surface-mapping) be used cautiously, for quick visualizations, and accompanied by careful confirmation that the mapping produced a reasonable result.

needed files

Before we can make a surface version of a volumetric image, we need to know what it's aligned to, so that we can pick the proper corresponding surface template. Recall that gifti surface files (*.surf.gii) are sort of the underlay anatomy for surface images (e.g., in my Getting Started post on Workbench we load surf.gii files to get a blank brain), so we'll need gifti surface files that will work with our volumetric image (and to serve as the underlay when we're ready to plot the converted volumetric ROI).

For this demo, we'll use fakeBrain.nii.gz (it should let you download without signing in; this is the same NIfTI as shown in other posts), which is aligned to MNI. One MNI dataset with the necessary surface files is the HCP 1200 Subjects Group Average Data; this post describes the files and gives the direct ConnectomeDB download link.

The HCP 1200 Subject Group Average download contains multiple *.surf.gii for each hemisphere, including midthickness, pial, and very_inflated. We can use any of these for visualization in Workbench, but which we pick for the volume to surface conversion does make a difference in what the resulting surface image will look like. It seems best to start with the midthickness surface for the conversion, then try others if the projection seems off.

using wb_command

The wb_command -volume-to-surface-mapping function does the conversion. wb_command.exe (on my Windows machine; file extension may vary) should be in the same directory as wb_view.exe, which you use to start the Workbench GUI. Don't double-click wb_command.exe - it's a command line program. Instead, open up a command prompt and navigate to the directory containing wb_command.exe (on my machine, /bin_windows64/). If you  type wb_command at the prompt it should print out some version and help information; if not, check if you're in the correct directory, and try ./wb_command if you're on a linux-type system.

Now we're ready: we give the function our input NIfTI (fakeBrain.nii.gz), our surface gifti (S1200.L.midthickness_MSMAll.32k_fs_LR.surf.gii), the output file we want it to make (demoL.shape.gii), and the options for it to use (-trilinear). Since surface gifti files are just for one hemisphere, we have to do the command twice, once for each. (I included the full path to each file below; update for your machine.)

 wb_command -volume-to-surface-mapping d:/temp/fakeBrain.nii.gz d:/Workbench/HCP_S1200_GroupAvg_v1/S1200.L.midthickness_MSMAll.32k_fs_LR.surf.gii d:/temp/demoL.shape.gii -trilinear  
 wb_command -volume-to-surface-mapping d:/temp/fakeBrain.nii.gz d:/Workbench/HCP_S1200_GroupAvg_v1/S1200.R.midthickness_MSMAll.32k_fs_LR.surf.gii d:/temp/demoR.shape.gii -trilinear  

We now have demoL.shape.gii and demoR.shape.gii, surface versions of fakeBrain.nii.gz, which can be viewed in Workbench (or other gifti-aware programs). Check the surface projection carefully: does the ROI align properly with the anatomy? If not, try a different .surf.gii or fitting option (e.g., -enclosing) in the wb_command call; these can make a big difference.

Below the jump are the surface images from the above commands, plotted on the S1200 midthickness in Workbench.

Holy crescents, Batman!

Quite a few of the posts over the last year or so have arisen from things that catch my eye as I review the SMS/MB4 images we're collecting in our ongoing project, and this is another. For quick comparison, I make (with knitr; we may give mriqc a try) files showing slices from mean, standard deviation, and tSNR images for participants, runs, and sessions.

Some participants have obvious bright crescent-shaped artifacts in their standard deviation images (the examples above are from two people; both calculated from non-censored frames, after completing the HCP Minimal Preprocessing pipeline). Looking over people and runs (some participants have completed 6 imaging sessions, over months), people have the crescents or not - their presence doesn't vary much with session (scanning day), task, or movement level (apparent or real).

They do, however, vary with encoding direction: appearing in PA phase encoding runs only. Plus, they seem to vary with subject head size, more likely in small-headed people (large-headed people seem more likely to have "ripples", but that's an artifact for another day).

All that (and thanks to conversations with practiCal fMRI and @DataLoreNeuro) gave a hint: these crescents appear to be N/2 ghost artifacts.

Playing with the contrast and looking outside the brain has convinced me that the crescents do align with the edges of ghost artifacts, which I tried to show above. These are from a raw image (the HCP Minimal Preprocessing pipelines mask the brain), so it's hard to see; I can share example NIfTIs if anyone is interested.

So, why do we have the bright ghosts, what should we do about it, and what does that mean for analysis of images we've already collected? Suggestions are welcome! For analysis of existing images, I suspect that these will hurt our signal quality a little: we want the task runs to be comparable, but they're not in people with the crescent: voxels within the crescent areas have quite different tSNR in the PA and AP runs.

Holy crescents, Batman! (We've been watching the 1966 Batman TV series.)

afni: 3dTstat with not-censored timepoints only

In a few recent posts I've shown images of the mean and standard deviation (calculated across time for each voxel), for QC tests. These are easy to calculate in afni (example here), but the 3dTstat command I used includes all timepoints (TRs), unless you specify otherwise. As described previously, we've been using a threshold of FD > 0.9 for censoring high-motion frames before doing GLMs. Thus, I wanted to calculate the mean and standard deviation images only including frames that were not marked for censoring (i.e., restrict the frames used by 3dTstat). This was a bit of a headache to code up, so R and afni code are after the jump, in the hopes it will be useful for others.