Summary of method

The statistical analysis of the fMRI data was based on a linear model with correlated errors. The design matrix of the linear model was first convolved with a hemodynamic response function modeled as a difference of two gamma functions (Glover, 1999) timed to coincide with the acquisition of each slice. Drift was removed by adding polynomial covariates in the frame times, up to degree 3(?), to the design matrix. The correlation structure was modeled as an autoregressive process of degree 1 (Bullmore et al., 1996). At each voxel, the autocorrelation parameter was estimated from the least squares residuals using the Yule-Walker equations, after a bias correction for correlations induced by the linear model. The autocorrelation parameter was first regularized by spatial smoothing with a 15(?)mm FWHM Gaussian filter, then used to `whiten' the data and the design matrix. The linear model was then re-estimated using least squares on the whitened data to produce estimates of effects and their standard errors.

In a second step, runs, sessions and subjects were combined using a mixed effects linear model for the effects (as data) with fixed effects standard deviations taken from the previous analysis. This was fitted using ReML implemented by the EM algorithm. A random effects analysis was performed by first estimating the the ratio of the random effects variance to the fixed effects variance, then regularzing this ratio by spatial smoothing with a 15(?)mm FWHM Gaussian filter. The variance of the effect was then estimated by the smoothed ratio multiplied by the fixed effects variance to achieve higher degrees of freedom.

The resulting T statistic images were thresholded using the minimum given by a Bonferroni correction and random field theory (Worsley et al, 1996).

References

Bullmore, E.T., Brammer, M.J., Williams, S.C.R., Rabe-Hesketh, S., Janot, N., David, A.S., Mellers, J.D.C., Howard, R. and Sham, P. (1996). Statistical methods of estimation and inference for functional MR image analysis. Magnetic Resonance in Medicine, 35:261-277.

Glover, G.H. (1999). Deconvolution of impulse response in event-related BOLD fMRI. NeuroImage, 9:416-429.

Worsley, K.J., Marrett, S., Neelin, P., Vandal, A.C., Friston, K.J. and Evans, A.C. (1996). A unified statistical approach for determining significant signals in images of cerebral activation. Human Brain Mapping, 4:58-73.

PowerPoint presentation of fMRIstat

The latest release of fmristat

The latest release has been written up in two papers submitted to NeuroImage: A General Statistical Analysis for fMRI Data and Estimating the delay of the hemodynamic response in fMRI data.

All programs now read and write ANALYZE format (if the extension is .img) as well as MINC format files (if the extension is .mnc). Fmrilm has been extended to estimate delays as well as effects with the addition of one more parameter, and fmridesign has been been extended to provide the necessary extra derivatives of the hrf. See the end of this document for instructions. (Note that the EXCLUDE parameter is no longer necessary for fmridesign and it must be ommitted.) The program can now estimate delays in seconds, shown below colouring the outer surface of the activated regions in top and cut-away front views. If you don't ask for delays, fmrilm should be completely backwards compatible with previous versions. Fmrilm has been merged with fmrilm_arp so that fmrilm can now fit autoregressive models of arbitrary order.

Multistat has been updated to fit a random effects model by REML using the EM algorithm. This gives a better analysis than the previous version for combining subjects (the previous version was better at combining runs on the same subject, but this is usually not as important as combining subjects).

Stat_threshold has been extended to find thresholds for F statistics as well as T statistics, and to find P-values for both if peak values are given instead. It also calculates P-values for the spatial extent of clusters of voxels above a threshold.

Locmax is a new function that finds local maxima of a minc file above a threshold.

Extract is a new function that extracts values from a minc file at specified voxels (in `register' format).

Example contains all the matlab commands for the example used here.

Data contains all the data for the example used here.

Plotting the hemodynamic response function (hrf) using fmridesign

1. PEAK1: time to the peak of the first gamma density;
2. FWHM1: approximate FWHM of the first gamma density;
3. PEAK2: time to the peak of the second gamma density;
4. FWHM2: approximate FWHM of the second gamma density;
5. DIP: coefficient of the second gamma density. The final hrf is: gamma1/max(gamma1)-DIP*gamma2/max(gamma2) scaled so that its total integral is 1.
If PEAK1=0 then there is no smoothing of that event type with the hrf. If PEAK1>0 but FWHM1=0 then the design is lagged by PEAK1. The default, chosen by Glover (1999) for an auditory stimulus, is:

Making the design matrices using fmridesign

1. EVENTID - an integer from 1:(number of events) to identify event type;
2. EVENTIMES - start of event, synchronized with frame and slice times;
3. DURATION (optional - default is 0) - duration of event;
4. HEIGHT (optional - default is 1) - height of response for event.

For each event type, the response is a box function starting at the event times, with the specified durations and heights, convolved with the hemodynamic response function (see above). If the duration is zero, the response is the hemodynamic response function whose integral is the specified height - useful for ‘instantaneous’ stimuli such as visual stimuli. The response is then subsampled at the appropriate frame and slice times to create a design matrix for each slice, whose columns correspond to the event id number. EVENT_TIMES=[] will ignore event times and just use the stimulus design matrix S (see next).

This is a sample run on Brian Ha's data, which was a block design of 3 scans rest, 3 scans hot stimulus, 3 scans rest, 3 scans warm stimulus, repeated 10 times (120 scans total). The hot event is identified by 1, and the warm event by 2. The events start at times 9,27,45,63... and each has a duration of 9s, with equal height:

Analysing one run with fmrilm

1: T statistic image, OUTPUT_FILE_BASE_mag_t.mnc or .img (for contrasts in the delays, mag is replaced by del in the extension). The degrees of freedom is DF. Note that tstat=effect/sdeffect. 2: effect image, OUTPUT_FILE_BASE_mag_ef.mnc or .img. 3: standard deviation of the effect, OUTPUT_FILE_BASE_mag_sd.mnc or .img. 4: F-statistic image, OUTPUT_FILE_BASE_mag_F.mnc or .img, for testing all rows of CONTRAST that deal with magnitudes (see CONTRAST_IS_DELAY below). The degrees of freedom are [DF, P]. 5: the temporal autocorrelation(s), OUTPUT_FILE_BASE_rho.mnc or .img. 6: the residuals from the model, OUTPUT_FILE_BASE_resid.mnc or .img, only for the non-excluded frames (warning: uses lots of memory). 7: the whitened residuals from the model, OUTPUT_FILE_BASE_wresid.mnc or .img, only for the non-excluded frames (warning: uses lots of memory). 8: the AR parameter(s) a_1 ... a_p, in OUTPUT_FILE_BASE_A.mnc or .img.

If WHICH_STATS to a row vector of length 8, then it is used for all contrasts; if the number of columns is less than 8 it is padded with zeros. Default is 1.

To get t statistics, effects and their sd's, and the F statistic for testing for any effect of either hot or warm vs baseline, and the autocorrelation:

To see the brain, use the first slice of the fMRI data > 500 as a mask:

To look at slices 4 to 6 simultaneously, try

F-tests should only be used when we are interested in any linear combination of the contrasts. For example, an F-test would be appropriate for detecting regions with high polynomial drift, since we would be interested in either a linear, quadratic or cubic trend, or any linear combination of these. The drift terms appear after the regressor terms, in the order constant, linear, quadratic, cubic etc. So to test for a drift, use the contrast

A linear effect of temperature

To model a quadratic effect of temperature (Figure b), add a third event type with a height equal to the temperature²:

Combining runs/sessions/subjects with multistat

For a more elaborate analysis, e.g. comparing the first two sessions with the next two, you can do it by:

Fixed and random effects

The ratio of the random effects standard error, divided by the fixed effects standard error, is outputed as ha_multi_sdratio.mnc. If there is no random effect, then sdratio should be roughly 1. If there is a random effect, sdratio will be greater than 1. To allow for the random effects, we multiply the sdeffect by sdratio, and divide the tstat by sdratio. The only problem is that sdratio is itself extremely noisy due to its low degrees of freedom, which results in a loss of sensitivity. One way around this is to assume that the underlying sdratio is fairly smooth, so it can be better estimated by smoothing sdratio. This will increase its degrees of freedom, thus increasing sensitivity. Note however that too much smoothing could introduce bias.

The amount of smoothing is controlled by the final parameter, fwhm_varatio. It is the fwhm in mm of the Gaussian filter used to smooth the ratio of the random effects variance divided by the fixed effects variance. 0 will do no smoothing, and give a purely random effects analysis. Inf will do complete smoothing to a global ratio of one, giving a purely fixed effects analysis (which we did above). The higher the fwhm_varatio, the higher the ultimate degrees of freedom of the tstat image (printed out before a pause), and the more sensitive the test. However too much smoothing will bias the results. The program prints and returns the df before doing any analysis; if it is too low, press control-c to cancel, and try increasing fwhm_varatio. The default is 15, and this is the recommended value. To invoke this, just run

Note that it is roughly 1 outside the brain and at most places inside, indicating little random effect over the runs, all the random effect sd is about 60% higher than the fixed effect sd (sdratio = 1.6) in posterior regions. The sdratio is much higher (~3) over subjects.

Thresholding the tstat image with stat_threshold and fdr_threshold

Fdr_threshold calculates the threshold for a t or F image for controlling the false discovery rate (FDR). The FDR is the expected proportion of false positives among the voxels above the threshold. This threshold is higher than that for controlling the expected proportion of false positives in the whole search volume, and usually lower than the Bonferroni threshold (printed out as BON_THRESH) which controls the probability that *any* voxels are above the threshold. The threshold depends on the data as well as the search region, but does not depend on the smoothness of the data. We must specify a search volume, here chosen as the first frame of the raw fMRI data thresholded at a value of 500:

References:

Locating peaks with locmax

Extracting values from a minc file using extract

To look at the time course of the data for the first run, and compare it to the fitted values, try this:

Estimating the time course of the response

The advantage of this is that it is much more flexible; it can be applied to event related designs, with randomly timed events, even with events almost overlaping, and the duration of the `events' need not be equal to the TR. Moreover, the results from separate runs/sessions/subjects can be combined using multistat.

Now all the time course is covered by events, and the baseline has disappeared, so we must contrast each event with the average of all the other events, as follows:

The modeled response and the fitted response are quite close. Note that the graph `wraps around' in time.

Estimating the delay

CONTRAST_IS_DELAY is a logical vector indicating if the row of CONTRAST is a contrast in the delays (1) or magnitudes (0). Note that you cannot estimate delays of the polynomial terms or confounds, and that if you estimate delays (1) then each non-zero element of CONTRAST must be estimable i.e. reference delay >0 for that response. If the length of CONTRAST_IS_DELAY is less then the number of contrasts, it is padded with zeros. Note that the t test for contrasts with a single number tests that the delay is equal to the reference delay (not that it equals zero). Default is 0, i.e. all contrasts refer to magnitudes Here is an example to estimate the delays of the hot, the warm and the difference hot-warm stimuli for the above data. In this case every contrast is a delay, so CONTRAST_IS_DELAY = [1 1 1]:

The first image is the delay, in seconds, where the t-statistic for the magnitude is greater than 4. The brain is outlined by thesholding the first fMRI scan at 500. The second image gives the sd of the delay. The delay is roughly 5.5 +/- 0.6 seconds, close to the nominal delay of 5.4 seconds for the Glover HRF.

Note that the tstat image is not quite the same as what you would get from fmrilm without asking for delays, i.e. with REFERENCE_DELAY=0. The reason is that in order to calculate delays, fmrilm fits a model with convolutions of the stimulus with the derivative of the hrf added to the linear model.

You could now put delay effects and sdeffects from separate runs/sessions/subjects through multistat as before.

Effective connectivity of all voxels with a reference voxel

Note that voxels in the neighbourhood of the pre-chosen voxel are obviously highly correlated due to smoothness of the image. To avoid very large t statistics, the t statistic image has a ceiling of 100 (the F statistic image has a ceiling of 1000).

Local maxima of this image could be identified, their values extracted and added as columns to ref_data to find further voxels that are correlated with the first reference voxel removing the effect of the second reference voxel. This could be repeated to map out the connectivity.

The same idea could be used with multistat to look for connectivity across runs or subjects; simply add an extra column to the design matrix X that contains the subject or run values at the reference voxel, and set the contrast to choose that column by puting a value of 1 for that column and zero for the others.

Higher order autoregressive models