FWHM & RESEL details for SPM and FSL

Executive Summary

Random Field Theory depends on a measure called the RESEL count to determine corrected P-values. How RESELs are reported on in SPM and FSL, however, differs. This post reviews the basics of RESELs and FWHM and details how to interpret the values reported by each software.

Background

Random Field Theory (RFT) is the method used in both FSL and SPM to make inference on statistic images. Specifically, it finds thresholds that control the chance of false positives while searching the brain for activations or group differences. The magic of RFT, unlike other methods, like Bonferroni, is that it accounts for the spatial smoothness of the data; as a result, a less severe correction is applied to highly smooth data while a more stringent correction is applied to relatively rough data.

The details of how RFT measures smoothness are technical (it is the inverse square root of the determinant of the variance-covariance matrix of the gradient of the component fields), but fortunately Keith Worsley introduced the notion of FWHM and RESELs [1] to make matters simplier.

For RFT, the FWHM is defined as the Full Width at Half Maximum of the Gaussian kernel required to smooth independent, white noise data to have the same “smoothness” as your data. Here, “smoothness” refers to the technical definition involving the “square root of the determinant…”. Crucially, FWHM is not the applied smoothness, e.g. the Gaussian kernel smoothing applied as part of preprocessing. It is the smoothness of the data fed into the GLM, which is a combination of the intrinsic smoothness of the data (affected by things like image reconstruction parameters and physiological noise). As we have 3-D data, FWHM is usually represented as a 3-vector, [FWHMx FWHMY FWHMZ], though be careful to check whether the units are in voxels or mm.

RESEL stands for RESolution ELement, and is a virtual voxel with dimensions [FWHMx FWHMY FWHMZ]. The RESEL count is the number of RESELs that fit into your search volume; in math mode:

R = V / ( FWHMX × FWHMY × FWHMZ )

where V is the search volume; note that V can be in voxels or mm3 (or whatever) as long as FWHMX, FWHMY & FWHMZ have the same units.

Note a possible source of confusion: A “RESEL” is 3D object, a cuboid; the “RESEL count” is a scalar, your search volume in units of RESELs.

SPM

SPM provides the most complete reporting of FWHM and RESEL-related quantities. In the Results figure, in the footer of the P-value table, SPM lists [FWHMx FWHMY FWHMZ] in both voxel and mm units. It also lists the search volume in RESELs (along with the same in mm3 and voxel units). Finally, it reports the size of one RESEL measured in voxels.

Note, if you check SPM’s arithmetic, you may find that R is not exactly V / ( FWHMX × FWHMY × FWHMZ ). The reason is that a more sophisticated method for measuring the search volume is used that involves counting voxel edges and faces, etc; see spm/src/spm_resels_vol.c & [2] for more details.

FSL

FSL provides no information about RESELs or FWHM in its HTML output. There is some information available, however, buried in the *.feat directories.

In each *.feat directory there is a stats subdirectory. In that directory you’ll find a file called smoothness that contains three lines, one each labeled DLH, VOLUME, RESELS.

DLH is the unintuitive parameter that describes image roughness; precisely it is

DLH = (4 log(2))3/2 / ( FWHMX × FWHMY × FWHMZ)

where FWHMX, FWHMY and FWHMZ are in voxel units.

VOLUME is the search volume in units of voxels. And, as a source of great possible confusion, RESELS is not the RESEL count, but rather, the size of one RESEL in voxel units

RESELS = FWHMX × FWHMY × FWHMZ

Sadly, the FWHMX, FWHMY and FWHMZ are never individually computed or saved. Note, however, the geometric mean of the FWHM’s is available as

AvgFWHM = RESELS1/3

which is some consolation.

Update: The verbose option to smoothest, -V, it will report the individual FWHMx, FWHMx, & FWHMy. But, again, this isn’t saved as part of feat output.

Who cares?

Any neuroimaging researcher should always check the estimated FWHM of their analysis. The reason is that if, for some reason the FWHM smoothness is less than 3 in any dimension, the accuracy of the RFT results can be very poor [3]. While the FSL user can’t examine the 3 individual FHWM, they can at least compute AvgFWHM and check this.

Another reason is to understand differences in corrected P-values. For example, say you conduct 2 studies each with 20 subjects. You notice that, despite having similar T-values, the two studies have quite different FWE corrected P-values. The difference must be due to different RESEL counts, which in turn can be explained by either a difference in search volume (V) or smoothness (FWHM).

Finally, one more reason to check the RESELs is when you are getting bizarre results, like insanely significant results for modest cluster sizes. For example, if VBM results are poorly masked it can happen that the analysis includes a large boundary of air voxels outside the brain. Not only is this unnecessarily increasing your search volume (possibly increasing your multiple testing correction), but the smoothness estimate in air will be dramatically lower than in brain tissue, and thus corrupt the accuracy of the inferences. A bizarre value for RESEL size or FWHM will reveal this problem.

Further Reading

For a light-touch introduction to Random Field Theory, see [3]; for more detail (but less than the mathematical papers) see the relevant chapters in these books [4,5,6].

References

[1] Worsley, K. J., Evans, A. C., Marrett, S., & Neelin, P. (1992). A three-dimensional statistical analysis for CBF activation studies in human brain. Journal of Cerebral Blood Flow and Metabolism, 12(6), 900–918. Preprint

[2] Worsley, K. J., Marrett, S., Neelin, P., Vandal, A. C., Friston, K. J., & Evans, A. C. (1996). A unified statistical approach for determining significant signals in images of cerebral activation. Human brain mapping, 4(1), 58–73.

[3] Nichols, T. E., & Hayasaka, S. (2003). Controlling the familywise error rate in functional neuroimaging: a comparative review. Statistical Methods in Medical Research, 12(5), 419–446.

[4] Poldrack, R. A., Mumford, J. A., & Nichols, T. E. (2011). Handbook of fMRI Data Analysis (p. 450). Cambridge University Press.

[5] Penny, W. D., Friston, K. J., Ashburner, J. T., Kiebel, S. J., & Nichols, T. E. (2006). Statistical Parametric Mapping: The Analysis of Functional Brain Images (p. 656). Academic Press. [Previous version available free]

[6] Jezzard, P., Matthews, P. M., & Smith, S. M. (2002). Functional MRI: An Introduction to Methods (p. 432). Oxford University Press.

Standardizing DVARS

Writing about web page http://www.nisox.org/Software/Scripts/FSL/#DVARS

Multiple authors have found the spatial standard deviation of the temporal difference data useful for detecting anomalous or artifactual scans in fMRI. Called DVARS by Power et al. (NeuroImage, 59(3), 2142-54), this measure is useful but lacks interpretable units. In particular, the nominal value of DVARS reflects both the temporal standard deviation and the degree of temporal autocorrelation.

I’ve created a simple approach to standardize DVARS, attempting to remove the dependence on temporal standard deviation and autocorrelation. See my notes on Standardizing DVARS and a FSL script DVARS.sh to implement this new version of DVARS.

SPM plot units

What are the units of a plot in SPM?

In SPM, when you click on “Plot” and ask for a “Contrast estimates and 90% CI”, what exactly are the units of the thing plotted? I get this question so many times I wanted to record my answer once and for all.

The short answer is: 2-3 times an approximation to percent BOLD change, if you have a long-block design; otherwise, it depends. Here’s the long answer.

The final units of a plotted contrast of parameter estimates depends on three things:

  1. The scaling of the data.
  2. The scaling of the predictor(s) that are involved in the selected contrast.
  3. The scaling of the selected contrast.

The scaling of the data

For a first level fMRI analysis in SPM the scaling is done “behind the scenes”… there are no options to change the default behavior: The global mean intensity for each image is computed, and the average of those values (the “grand grand mean”) is scaled to be equal to be 100. However, SPM’s spm_global produces a rather dysmal estimate of typical intracerebral intensity (note that the author is listed as “Anonymous”; to be fair, it worked well on the tightly-cropped, standard-space PET images that it was originally designed for). The estimate of spm_global is generally too low, and so when the data are divided by it and multiplied by 100, the typical brain intensity is well over 100. In my experience, if you examine a beta_XXXX.img corresponding to a constant in a first-level regression model, typically gray matter values range from 200 to 300 when it should be 100; for first-level analyses done in subject-space (where there are even more air voxels) the bias can be even greater. Let’s call the typical gray matter intensity of your BOLD T2* images to be B for later correcting this effect (e.g. B = 250 or 300 or whatever you find).

At the second level, you depend on the contrast estimates passed up from the first level to have a valid interpretation. So all three of these rules should apply to the first and second level models, but the data scaling is done at the first level (no further scaling should be done at the second level; e.g., in the Batch Editor, second-level models have “Global normalisation -> Overall grand mean scaling” set to “No”).

The scaling of the predictor(s) that are involved in the selected contrast

At the first level, the predictors must be scaled appropriately to give the regression coefficients the same units as the data. For BOLD fMRI, this means the baseline to peak magnitude should be 1.0… i.e. a one unit change in a beta coefficient will result in a one unit change in predicted BOLD effect. In SPM, for loooong blocks/epochs, 20 seconds or longer, the scaling is as expected, with the baseline-to-plateau difference being 1.0. For shorter blocks however, the peak can be higher – as this corresponds to the ‘initial peak’ part of the HRF, as the duration is too short for the flat plateau to be established (this ‘early peak’ is due to symmetry of convolution, as it must exactly mirror the post-stimulus overshoot; there is some physiological basis for this as well, though) – or they can be lower. For events, the HRF is normalized to sum to unity, which doesn’t help you predict the peak.

If you have anything but long blocks, then, you have to examine the baseline-to-peak (or if you prefer, baseline-to-trough) height, and make a note of it for later scaling (call it H). (The design matrix is in SPM.xX.X, in SPM.mat; plot individual columns to judge the scaling). For event-related designs, however, it can be tricky to uniquely define the baseline-to-peak height; if you have a slow design, with isolated events, you’ll find a different peak height than if you have a design with some closely-spaced events as the responses pile-up additively. A reasonable unique definition of a event-related response magnitude is the height of an isolated event; again, call this H.

At the second level, you must also make sure that your regressors preserve the units of the data. For example, if you are manually coding dummy regressors, like gender, as -1/1, realize that the unit effect of beta coefficient for gender will result in a two unit change in the data. Safer is using a 0/1 coding, though that then re-defines the intercept; best yet is to use -1/2 and 1/2 as dummy values for a dichotomous variable.

For parametric modulation and covariates, note that a regression coefficient’s units are determined by the units of the covariate. For example, if you have a parametric modulation in a reward experiment, where the units are £’s that are at risk in a bet, then your regressor coefficient has units (if everything else is working as it should) “percent BOLD change per £ at risk”. Likewise, if you use an age predictor in the second level, if you enter age in years (centered or uncentered), the units of the age coefficient are “percent BOLD change per year”

Finally, beware of trying to compare between percent BOLD change in Block Designs vs. Event Related designs. Due to the nonlinearities of the BOLD effect, a percent BOLD change for an event related response is really a different beast from a percent BOLD change for an block design. When reporting one type, be sure it is clear what sort of effect you’re reporting, and don’t try averaging or comparing their magnitudes.

The scaling of the selected contrast

The General Linear Model is a wondrous thing, as you’ll get the same T-values and P-values for the contrast [-1 -1 1 1] as you do for [-1/2 -1/2 1/2 1/2]. However, the interpretation of the contrast of parameter estimates (contrast times beta-hats) is different for each. The following rules to scale contrasts will usually work to ensure that a contrast preserves the units of the regression coefficients:

  • Sum of any positive contrast elements is 1.
  • Sum of any negative contrast elements is -1.

If the contrast elements are all positive or all negative, the interpretation is clear, as you are ensuring that the beta-coefficents (each which should have unit-change interpretations) will be averaged, giving a (average) unit change interpretation. If there are a mix of positive and negative elements, I see this rule as carefully ensuring there is a positive unit change quantity subtracted from a negative unit change quantity, giving the expected difference of unit change. (This rule should work for covariates and simple 0/1 dummy variables; this FSL example shows a case with -1/0/1 dummy variables where the contrasts violate these rules yet actually give the correct units.)

When it goes right, when it goes wrong

If everything had gone well, then the plotted units are “approximate percent BOLD change”. The “approximate” qualification results from only using a global-brain mean of 100; an “exact” percent BOLD change requires that the voxel in question has had baseline (or grand mean) intensity of exactly 100. We don’t grand-mean-normalise voxel-by-voxel because of edge effects: At the edge of the brain, the image intensity falls off to zero, and dividing by these values will give artifactually large values.

Of the three scalings in play, the first two of these three items are basically out of the control of the user. As hinted, first-level fMRI data in SPM are generally mis-scaled, resulting in global brain intensity between 200 and 300, i.e., 2 to 3 times greater than anticipated. Also, if you don’t have long blocks, you’ll need to apply a correction for the non-unity H you measured. The corrections are as follows: Dividing the reported results B/100 (i.e. multiplying by 100/B) will correct the problem of the data being incorrectly scaled. Multiplying by H will correct the units of the regression coefficients. After both of these corrections, you should get something closer to approximate percent BOLD change.

This all sounds rather depressing, doesn’t it? Well, there is a better way. If you use the Marsbar toolbox it will carefully estimate baseline intensity for the ROI (or single-voxel ROI—make sure it isn’t on the edge of the brain though!) and gives exact percent BOLD change units.

Postscript: What about FSL?

FSL has an good estimate of the global mean, because it uses morphological operations to identify intracerebral voxels. It scales grand mean brain intensity to 10,000, however (because it used to use a integer data representation and 100 would have resulted in too few significant digits.) Block-design predictors are scaled to have unit height; the default HRF doesn’t use a post-stimulus undershoot, and so baseline-to-peak height is the same as trough-to-peak height. Event-related predictors are not scaled to have peak height of unity; however, the recommended tool for extracting percent BOLD change values, featquery, estimates a trough-to-peak height and then appropriate scales the data. However, as Jeanette Mumford has pointed out, this isn’t really the right thing to do when you have closely spaced trials; see her excellent A Guide to Calculating Percent Change with Featquery to see how to deal with this.

Other references

Thanks to Jeanette Mumford & Tor Wager for helpful input and additions to this entry, and Guillaume Flandin for the additional references.

Random Shell Tricks

This is quite boring… just a list of useful shell commands that I’m sure I’ll need again sometime.

  • Find files/directories in a specified date range, do something to them; base this on the top level only, not traversing into each subdirectory.

    find ./* -prune -newermt 1/1/2010 -not -newermt 1/1/2011 -exec mv {} ../2010/ \;

    will move all files/directories with modification times in 2010 to a directory one level up, called 2010. I use this to clean up a download/scratch/temp directory (I should just delete the files, but just in case, I organize them into yearly folders).

FSL Scripts

Writing about web page http://www.nisox.org/Software/Scripts/FSL

I have recently uploaded a number of FSL scripts which folks might find useful. See my FSL scripts page for details, but the new ones include

  • fslstats.sh Allows you to use fslstats with multiple files
  • fsl_fdr.sh Gives a much nicer interface to FDR results in FSL
  • Dummy.sh Creates dummy variables, handy for manual/scripted creation of design matrices.

Creating Paired Differences

If you analyze 2-time point longitudinal data, you will eventually observe that it is easier to create and analyze difference data (e.g. data from time 2 minus data from time 1). However, if you’re not a scripting guru this might be an annoying task.

This PairDiff.m script will take an even number of images and produce a set of paired difference images, even minus odd.

function PairDiff(Imgs,BaseNm)
% PairDiff(Imgs,BaseNm)
% Imgs   - Matrix of (even-number of) filenames
% BaseNm - Basename for difference images
%
% Create pairwise difference for a set of images (img2-img1, img4-img3, etc),
% named according to BaseNm (BaseNm_0001, BaseNm_0002, BaseNm_0003, etc).
%______________________________________________________________________
% $Id: PairDiff.m,v 1.2 2012/02/16 21:36:41 nichols Exp $

if nargin<1 | isempty(Imgs)
  Imgs = spm_select(Inf,'image','Select n*2 images');
end
if nargin<2 | isempty(BaseNm)
  BaseNm = spm_input('Enter difference basename','+1','s','Diff');
end
V = spm_vol(Imgs);
n = length(V);
if rem(n,2)~=0
  error('Must specify an even number of images')
end

V1 = V(1);
V1  = rmfield(V1,{'fname','descrip','n','private'});
V1.dt = [spm_type('float32') spm_platform('bigend')];

for i=1:n/2
  V1.fname = sprintf('%s_%04d.img',BaseNm,i);
  V1.descrip = sprintf('Difference %d',i);
  Vo(i) = spm_create_vol(V1);
end

%
% Do the differencing
%

fprintf('Creating differenced data ')

for i=1:2:n

  img1 = spm_read_vols(V(i));
  img2 = spm_read_vols(V(i+1));

  img = img2-img1;

  spm_write_vol(Vo((i+1)/2),img);

end

After I created this I realized that Ged Ridgway also has a similar script, make_diffs.m , that takes two lists of images (baseline, followup) and does the same thing, though with perhaps more intuitive filenames.

Flame without 1st level directories – copes only.

Follow-up to Flame without 1st level directories from the NISOx blog (formerly Neuroimaging Statistics Tips & Tools)

My earlier post described how to get the basic FEAT second-level results when you just have a stack of COPE and VARCOPE images. In fact, you don’t even need the VARCOPE images.

If you just have the COPE images, in the flameo call omit --vc=4dvarcope (which you don’t have) and replace the --runmode=flame1 with --runmode=ols and you should be all set.

Save As HTML in Powerpoint 2011 for Mac

Follow-up to Extracting Images from Powerpoint Presentations for (e.g.) LaTeX presentations from the NISOx blog (formerly Neuroimaging Statistics Tips & Tools)

[NOTE: Powerpoint can now do this directly via “Export as…”, and then selecting HTML. Leaving this entry just for reference. -TN 2017/04/27]

Arrrgh. In Powerpoint 2011 for Mac there is no “HTML” option in the “Save as…” dialog.

Fortunately, there is a Microsoft Support note that gives a work around. For convenience (& my memory) here it is distilled:

  1. Press ALT+F11. (Depending on your keyboard settings, you might have to hold down “fn” to get F11.) This will open a “Project” window.
  2. Press CMD+CTRL+G to open the “Immediate” window. (The support page leaves out CMD, the ‘Apple’ key.)
  3. In the “Immediate” window paste in the following:
    ActivePresentation.SaveAs "/tmp/name.htm", ppSaveAsHTML, msoFalse
    Press ENTER.

This will give you a name.htm file and a name_files directory. Of course, you can replace ”/tmp/myfile.htm” with the desired path and filename. However, I find it saves keystrokes to save it into /tmp, and then move the file to wherever you really want it.

While you’ll see the images in the *_files/ folder, disappointingly, I don’t find the same original source images like I have seen before. Again, Arrrgh.

-Tom

SPM8 Gem 1: Zero NaN’s with the zeronan.m script

Follow-up to SPM99 Gem 3: NaNing zero values from the NISOx blog (formerly Neuroimaging Statistics Tips & Tools)

This was the topic of SPM99 Gem 3, converting NaN’s to zeros. For SPM8, see the following script zeronan.m that will zero NaN’s for you.

-Tom

function ofNm = zeronan(ifNm,val)
% FORMAT ofNm = zeronan(ifNm,val)
% ifNm  - Input filename(s)
% val   - Value to set NaN's to (defaults to zero)
%
% Output:
% ofNm  - Cell array of output filenames.
%
%
% Images have NaN's replaced with zeros, and new images, prefixed with a
% 'z', are created.
%
%________________________________________________________________________
% Based on zeronan.m,v 1.3 2005/10/26 21:58:55 nichols Exp 
% Thomas Nichols, 1 April 2011

if nargin<2, val = 0; end
if nargin0'); end

if ~iscell(ifNm)
  ifNm = cellstr(ifNm)';
else
  ifNm = ifNm(:)';
end

OtfNm = {};

for fNm = ifNm

  fNm = fNm{:};

  OfNm = ['z' fNm];
  [pth,nm,xt,vol] = spm_fileparts(fNm);
  OfNm = fullfile(pth,['z' nm xt]);

  % Code snippet from John Ashburner...
  VI       = spm_vol(fNm);
  VO       = VI;
  VO.fname = OfNm;
  VO       = spm_create_vol(VO);
  for i=1:VI.dim(3),
    img      = spm_slice_vol(VI,spm_matrix([0 0 i]),VI.dim(1:2),0);
    tmp      = find(isnan(img));
    img(tmp) = val;
    VO       = spm_write_plane(VO,img,i);
  end;

  OtfNm = {OtfNm{:} OfNm};

end

if nargout>0
  ofNm = OtfNm;
end

Extracting Images from Powerpoint Presentations for (e.g.) LaTeX presentations

Have you ever needed to create a LaTeX presentation but needed figures that only lived in a PowerPoint file? Or have you found that the only copy of a crucial figure you need lives inside a PowerPoint presentation, and when you copy it out of PowerPoint into another application it looks horrific?

This is a cool trick to extract the original (or close to original) image file from the presentation.

Open the presentation, delete all the slides you don’t need. The option is different with different versions of Powerpoint; look for “File -> Export as…” or ” “File -> Save as Web Page…”. You will get a .html (or .htm) file, and a _files directory. Inside the _files directory you will find a bunch of image files. It might take some poking around if there are lots of images on the slide, but you should find a close to original version of the image pasted into the slide.