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 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. 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 regularizing 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.

References

PowerPoint presentation of fMRIstat

The latest release of fmristat

25 April, 2002: Fmrilm has been extended to allow for different confounds for different slices. This is useful for effective connectivity, where a new function fmri_interp has been added to resample a single fmri time series at different frame times and slice times. See Effective connectivity of all voxels with a reference.

12 April, 2002: Fmristat now reads ANALYZE format files created on any machine, and the file seperator can be / or \.

8 April, 2002: Fmrilm now reads multiple 3D image files as well as a single 4D file - useful for ANALYZE format data.

Fmristat now allows for non-isotropic (non-constant FWHM) data. Many of the FWHM parameters can now be replaced by FWHM images. Fmristat now works for 2D data (one slice) as well as 3D data (more than one slice).

Fmrilm and multistat have been extended to estimate the effective FWHM of the data both from fixed and random effects. This allows for non-isotropic inference.

Multistat has been extended to estimate mixed effects conjunctions. It is now very fast if no sdeffect files are specified, making it suitable for analyzing PET data.

Locmax now finds clusters as well as local maxima.

Stat_threshold has been extended to find P-values and thresholds for conjunctions, and it now adjusts P-values and thresholds for clusters if the cluster extent is measured in resels and the resels are themselves estimated by fmrilm or multistat .

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

Data contains all the data for the example used here.

New functions

Pca_image does an exploratory Principal Components Analysis (PCA).

Glass_brain produces a maximum intensity projection or 'glass brain'.

Fdr_threshold finds Q-values and thresholds that control the false discovery rate.

Stat_summary produces an SPM-stlye summary of clusters, peaks and their P-values and Q-values.

Conf_region produces approximate confidence regions for the location of peaks in a statistic image.

Gauss_blur blurs 3D images with a Gaussian filter.

Looking at the fMRI data using pca_image

The first component, with an sd of 105.7 that accounts for 77.8% of the variance of the data, has extreme temporal values in the first few frames, perhaps because the EPI scanner has not reached a steady state. The spatial component is roughly the difference between the first few frames and the average of the others. The second component appears to capture a linear drift. We can see this more clearly by excluding the first two frames as follows:

Now the first component appears to have captured the linear drift around the edges of the brain, with an sd of 27.8 that explains 22.9% of the variability of the fMRI data. Note that signs are arbitrary, e.g. we can change the signs of both the time and space components. This component is saved in the first column of V; we could use it later on as the CONFOUND parameter to fmrilm to remove drift, instead of, or in addition to removing polynomial drift.

The next two components, with sd's of 17.6 and 15.7, appear to capture a cyclic pattern with 10 cycles that might be related to the 10 repetitions of the heat stimulus paradigm. This cyclic pattern is expressed in the supplementary motor areas which is confirmed below by a model-driven analysis. The fourth component, with an sd of 12.1, appears to be noise.

Pca_image can also do a PCA on the coefficients from a linear model specified by a design matrix X_interest, after removing effects from another linear model specified by a design matrix X_remove (Worsley et al., 1998). As an example, we might try to confine the PCA to the 12 scans in each block. To do this, we set X_interest to a design matrix for the 12 time points within each block. We might also remove a constant and the linear drift already detected as the first component of the PCA above:

The stimulus was 3 frames each of rest, hot, rest, warm, repeated 10 times, but delayed 6s (2 frames) by the hemodynamic response. The first component appears to have captured a warm minus hot stimulus effect. The second component, though much smaller, appears to be capturing a 6s delay of the first component. The rest appear to be noise.

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.
9: the estimated FWHM in OUTPUT_FILE_BASE_FWHM.mnc or .img - note this uses more time and memory.

If WHICH_STATS to a row vector of length 9, then it is used for all contrasts; if the number of columns is less than 9 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, the autocorrelation, and the FWHM:

Visualizing the results using glass_brain

To see the brain using the first slice of the fMRI data as a mask:

To get an SPM-style maximum intensity projection ('glass brain') above 3 try:

F-tests

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

These two parameters are only used to calculate the final degrees of freedom, and don't affect the images very much. Finally, for a fixed effects analysis (not recommended!) we set the last parameter to Inf (see later). The final run is:

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

Converting the analysis to percent of raw fmri data using fmri2pcnt

Fmri2pcnt will also accept a matrix of file names to be converted. It will then convert each file to percent, and return a matrix of the names of the converted files. So to redo the analysis on percentages using multistat, try this:

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 efficiency. 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 efficiency. 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 default is 15, but a value of 20 mm gives a final df that is reasonably high. 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. 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.

The effective FWHM over runs looks very similar to that over scans, but a lot more noisy, so some smoothing with gauss_blur makes it easier to interpret:

Reference:

Thresholding the tstat image with stat_threshold and fdr_threshold

Often the FWHM is estimated as an image rather than a fixed value. The extra variability in this estimate increases P-values and thresholds. It does not have much effect on P-values and thresholds for local maxima if the search region is large, but it can have a strong effect on P-values and thresholds for spatial extent, since clusters are small. The random noise in the FWHM is measured by the degrees of freedom of the linear model used to estimate it. For example, if the FWHM comes from a single run of fmrilm, then the df of the estimated FWHM is 112. To incorporate this into stat_threshold, add an extra row to the df:

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 mask_thresh=452:

References:

Locating peaks and clusters with locmax

Producing an SPM style summary with stat_summary

Confidence regions for the spatial location of local maxima using conf_region

Reference:

Conjunctions

Note that the peaks and thresholds are a lot lower than for a conventional T statistic image; the P=0.05 threshold for a 1000000mm^3 brain as above is

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

The first image is the delay, in seconds, where the t-statistic for the magnitude is greater than 4.86. The brain is outlined by thesholding the first fMRI scan at mask_thresh=452. The second image gives the sd of the delay. The delay shify is roughly 0.5 +/- 0.6 seconds, indicating no significant delay shift from the Glover HRF.

If a contrast involves the delays of more than one response, e.g. the third contrast which compares the delays of the hot and warm response, then it only makes sense to look at places where the magnitudes of all the responses involved are large. To do this, try:

Note that if you ask for magnitudes as well as delays in the same analysis, i.e.

The best way of visualizing delays is to use register like this:

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

Efficiency and choosing the best design

You can also get the actual linear combination of the data that is used to estimate these effects. These weights show how the frames are actually contributing to the effect. The array Y (second component of the output) contains the weights for each frame (first dimension), contrast (second dimension) and frame (third dimension). For the first 40 frames of slice 4:

The hot weights don't quite match the hot response because they are adjusted for the warm response, which is not quite orthogonal (if they were orthogonal, the weights would match the responses). The weights of the hot-warm are in fact the hot weights minus the warm weights. The weights are also adjusted for the drift as well; this is more evident if we have longer blocks, e.g. the extreme case of two 90 second blocks instead of twenty 9 second blocks:

Now the sd of the effects are much larger (this is a less efficient design), particularly for the warm stimulus. Part of the reason is that the responses are now correlated (confounded) with the cubic drift. This can be seen more clearly in the plot of the weights, which are concentrated at the beginning and end of the blocks (particularly for the warm stimulus). In other words, because the long blocks look a bit like drift, most of the information about the response comes from the frames when the stimulus is changing rapidly. The steady state of the response in the middle of the blocks is not providing much information because it is indistinguishable from (condfounded with) drift. This also explains why this design is less efficient: the extra frames in the middle of the blocks are wasted because they don't contribute to the effect. This is reflected in the increased sd of the effects, as already noted.

efficiency can also be used to find the sd of delays as well as magnitudes. Interestingly, the sd is in seconds irrespective of the data. This sounds a little strange, but if you think about it, delay estimation does not depend on the scaling of the data. To do this, just set contrast_is_delay to 1 where you want sd's of delays instead of magnitudes. For a temporal correlation of 0.3 (typical of cortex), cubic drift, and no confounds:

By looping through a variety of choices of stimulus duration and interstimulus interval, you can find out which combination gives the best design. The code is given in example.m, but the result is:

(a) A stimulus duration and interstimulus interval of about 8 seconds seems to be optimal for estimating the hot stimulus (the same is true of the warm stimulus). Fortunatley any values in the range 6 to 16 seem to be close to optimal.
(b) A stimulus duration of about 9 seconds and NO interstimulus interval seems to be optimal for hot-warm - this makes sense, because there is no need for the rest periods between the hot and warm stimuli. Again anything in the range 5 to 20 seems to be equally good, and short stimuli should be avoided.
(c) For delays, shorter bloocks are optimal, attaining sd's of about 0.2 seconds. Note that for very short blocks the sd of the magnitude is too large to reliably estimate the delay (black areas).
(d) Again shorter blocks seem to be better for delay differences.

Event related designs can be handled by setting the duration to zero. For a fixed number of events, we tried three different designs: equally spaced events ('uniform'), randomly spaced events ('random'), and all the events concentrated together at a single time point in the middle of the frames ('concentrated'). This last design is an extreme case, equivalent to a single event with magnitude = number of events. We then varied the number of events, and plotted the sd of the effect against the average time between events = time interval for the experiment (360s) divided by the number of events. Again the code is in example.m and the results are shown here:

For the magnitude (solid lines), and well seperated events, uniform and random designs give much the same increase in sd as the time between events increases. This is because as the time between events increases, the number of events decreases, and the sd is proportional to 1/sqrt(number of events) if the events are well separated. The uniform design seems to be optimumal at about 12s between events and gets less efficient as the time between events gets shorter than 12s. This is because after blurring by the HRF, the 'rest' periods get shorter and so there is less baseline to compare the stimulus with. On the other hand, the random design continues to show improvement as the time between events decreases because there is always some baseline present even after blurring with the HRF.

For the delays (dottted lines), it is interesting to note that the uniform design appears to be better, whereas the random design is better for the magnitude. Presumably the more consistent spacing between uniform events is better for estimating the shape of the HRF, including the delay. Also shorter times between events are better, but the results for very short times may be unreliable because the response of closely spaced events may be less than that of widely spaced responses, i.e. the additivity of responses may break down. For example, 10 presentations of a visual stimulus 0.1s apart might evoke less total response than 10 presentations 10s apart. The same thing might happen for blocks: 1s of continuous pain might evoke less total response than 10 well-separated presentations of 0.1s of pain.

This is exagerated when we force all the events very close together at a single time point (red lines). The sd of the effect is now inversely proportional to the number of events concentrated at the point, which is itself inversley propotional to the average time between events, so the sd of the effect is directly proportional to the average time between events. According to the graph, we reach the absurd conclusion that this is by far the best strategy! This shows that the additivity assumption must be wrong.

Conclusion: For block designs, try 8 or 9s blocks; for event-related designs, try about 12s or more between events, either uniformly or randomly spaced. For delays, shorter blocks and times between events are better. Fortunately halving or doubling these values produces only slightly less efficient designs, so there is great flexibility in the design parameters.

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 10000).

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

The images of the two autoregressive coeficients show that the second one, a_2, is near to zero so the first order autoregressive model is adequate, which is what we were using before (NUMLAGS=1, the default).