Summary of method

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 Gaussian filter. The variance of the effect was then estimated by the smoothed ratio multiplied by the fixed effects variance. The amount of smoothing was chosen to achieve 100 effective degrees of freedom.

The resulting T statistic images were thresholded using the minimum given by a Bonferroni correction and random field theory, taking into account the non-isotropic spatial correlation of the errors.

References

PowerPoint presentation of fMRIstat

The latest release of fmristat

Looking at the fMRI data using pca_image

The first component, with an sd of 0.68 that accounts for 46.9% of the correlations 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, revealing a sharp drop in ventricles and cortex. The second component appears to capture a linear drift. We can see this more clearly by excluding the first three frames as follows (note that you don't need to specify mask_thresh, pca_image will call fmri_mask_thresh to find it):

Now the first component appears to have captured the linear drift around the edges of the brain, with an sd of 0.38 that explains 14.7% of the correlations 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 spline drift.

The next two components, with sd's of 0.28 and 0.24, 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 Partial Least Squares PLS) or model-driven analysis. The fourth component, with an sd of 0.19, appears to be noise or perhaps the 'resting state network'.

Pca_image can also do Partial Least Squares (PLS), i.e. a PCA on the (normalised) 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 do a PLS on the time course within each block by confining 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 late peaking hot minus warm stimulus effect very similar to the fmrilm analysis (see later). The second component appears to be capturing an earlier onset of the heat stimulus. 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 the test 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

_mag_t T statistic image =ef/sd for magnitudes. If T > 100, T = 100.
_del_t T statistic image =ef/sd for delays. Delays are shifts of the time origin of the HRF, measured in seconds. Note that you cannot estimate delays of the trends or confounds.
_mag_ef effect (b) image for magnitudes.
_del_ef effect (b) image for delays.
_mag_sd standard deviation of the effect for magnitudes.
_del_sd standard deviation of the effect for delays.
_mag_F F-statistic for test of magnitudes of all rows of CONTRAST selected by _mag_F. The degrees of freedom are DF.F. If F > 1000, F = 1000. F statistics are not yet available for delays.
_cor the temporal autocorrelation(s).
_resid the residuals from the model, only for non-excluded frames.
_wresid the whitened residuals from the model normalized by dividing by their root sum of squares, only for non-excluded frames.
_AR the AR parameter(s) a_1 ... a_p.
_fwhm FWHM information:
Frame 1: effective FWHM in mm of the whitened residuals, as if they were white noise smoothed with a Gaussian filter whose fwhm was FWHM. FWHM is unbiased so that if it is smoothed spatially then it remains unbiased. If FWHM > 50, FWHM = 50. Frame 2: resels per voxel, again unbiased. Frames 3,4,5: correlation of adjacent resids in x,y,z directions.

e.g. WHICH_STATS='try this: _mag_t _del_ef _del_sd _cor_fwhm blablabla' will output t for magnitudes, ef and sd for delays, cor and fwhm. You can still use 1 and 0's for backwards compatibility with previous versions - see help from previous versions.

Visualizing the results using view_slices, glass_brain and blob_brain

To see the brain using the first frame of the fMRI data as a mask (note view_slices slices start at 0, so subtract 1):

or to see the first 12 slices:

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

To get a 3D 'blob' plot of voxels where the T statistic > 5, coloured by the hot - warm effect (in %BOLD), try blob_brain:

Also of interest is the lag 1 autocorrelation smoothed by fwhm_cor=12.4444mm (see above), which is much higher (~0.3) in cortical regions

The estimated FWHM of the data is close to 6 mm outside the brain (due to motion correction), but much higher (~12mm) in cortical regions, averaging to fwhm_data=8.7855 in the whole brain (see above):

F-tests

Note that we can use stat_threshold to find the threshold for significant drift (P<0.05), which is ~11.39 (the F statistic has df.F = 3 95) (see later). Another good use of the F-test is for detecting effects when the hemodynamic response is modelled by a set of basis functions (see Estimating the time course of the response ).

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

Fixed and random effects

The ratio of the random effects standard error, divided by the fixed effects standard error, is outputted as test_multi_hmw_rfx.mnc. If there is no random effect, then rfx should be roughly 1. If there is a random effect, rfx will be greater than 1. The only problem is that rfx 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 rfx is fairly smooth, so it can be better estimated by smoothing rfx. 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, and the more sensitive the test. However too much smoothing will bias the results.

By giving fwhm_varatio a negative value, it treats this as a target df, and adjusts the smoothing to try to achieve this. The default is -100, i.e. multistat finds the amount of smoothing to achieve 100 df. Why 100? The idea is that even if the assumptions underlying the df formula are wrong by a factor of 2, this will not greatly affect the distribution of T statistics, since T statistics with more than 50 df are very close to Gaussian. Using this default:

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 30% higher than the fixed effect sd (rfx = 1.3) in posterior regions.

For combining over subjects, the rfx is much higher (~3) and less smoothing should be used otherwise rfx is biased and you get more false positives. We recommend less than 10mm smoothing (fwhm_varatio=10), which can usually be achieved by aiming for less df, e.g. 40 df (fwhm_varatio = -40).

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

Another way of detecting significant activation is by looking at the spatial extent of activated regions (clusters), rather than peak height. We first set a lower threshold, usually in the range 2.5 - 4. Lower thresholds are better at detecting large clusters, but setting it too low invalidates the P-values. By default, stat_threshold uses the threshold for P=0.001 (uncorrected) which is cluster_threshold = 3.17. The resulting threshold for the spatial extent of clusters is extent_threshold = 638 mm^3. This means that any cluster of neighbouring voxels above a threshold of 3.17 with a volume larger than 638 mm^3 is significant at P<0.05.

Sometimes we are interested in a specific peak or cluster e.g. the nearest peak or cluster to a pre-chosen voxel or ROI. We don't have to search over all peaks and clusters, so the thresholds are lower. The threshold for peak height is then peak_threshold_1 = 4.30, and spatial extent_threshold_1 = 221 mm^3 is the extent of a single cluster chosen in advance - see Friston (1997).

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. 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 equals the df of T statistic image. Suppose this is 100, for comparison (actually it was 99). 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:

Finding the exact resels of a search region with mask_resels

The *real* perfectionist might want to allow for non-isotropic random fields as well, i.e. non-uniform FWHM. This requires quite a bit more work. We first need to save normalised residuals from fmrilm or multistat (depending on which program generated the test statistic) then replace the FWHM (before 8mm) by a file name fwhm_info that contains these residuals:

Locating peaks and clusters with locmax

Producing an SPM style summary with stat_summary

A file with _cluster before the extension is also produced which takes as its values the cluster id, and 0 elsewhere. If you put this into e.g. 'register' you can see where the clusters are. A glass brain, together with a 3D rendering of the significant clusters (which you can rotate), is produced as well:

Finally, a file with _Pval before the extension is also produced which takes the negative of the corrected P-values (so that bigger values are more significant and it looks better in e.g. `register') with values of -1.1 outside the mask:

Confidence regions for the spatial location of local maxima using conf_region

Reference:

Conjunctions

We first form the combined images of the hot and warm stimuli as above, and producing the p-value images for each:

Now we find the conjunction (minimum) of the negative p-values using conjunction, which also produces a list of clusters and local maxima like stat_summary:

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

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 delays and their sd (second and third figures) only make sense where the magnitude of the response is significantly large (first figure). The delay shift is roughly 0.5 +/- 0.6 seconds, indicating no significant delay shift from the Glover HRF (fourth figure). You could also put delay effects and sdeffects from separate runs/sessions/subjects through multistat as before.

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. Unfortunately there are no points where both the hot and the warm stimulus are exceed the 4.93 threshold, so we can't compare the delays.

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:

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 3 (zero based):

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.

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

These analyses address the question of what voxels are correlated with the reference voxel, after subtracting out the stimulus, in other words, the so-called 'resting state network'. This may not be very interesting. The more intersting question is how the connectivity might be modulated by the pain stimulus, compared to the warm stimulus, that is, how the stimulus changes the connectivity. To do this, we must add an interaction (product) between the hot and warm stimuli, and the reference voxel time course. The function XC produces a new set of confounds, which add to the confounds the product of every event type in X_cache with every covariate in the confounds. Confound covariates run fastest, so the order in X_confounds is [cov1, cov2,... ,event1*cov1, event1*cov2,... event2*cov1, event2*cov2, ...]. Once again we can pick out the hot-warm contrast with contrast.C=[0 1 -1], i.e. 0 for the single covariate (reference voxel time course), then 1 and -1 to pick out the product of the hot minus warm stimulus with the covariate.

Unfortunately there is no evidence that the connectivity is modulated by the pain in this example.

Higher order autoregressive models

The images of the autoregressive coefficients show that the last three are near to zero so the first order autoregressive model is adequate, which is what we were using before (NUMLAGS=1, the default).