AFNI example

Looking at the fMRI data using pca_image

The first component 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 4 frames as follows (note that you don't need to specify mask_thresh, pca_image will call fmri_mask_thresh to find it):

We might try to confine the PCA to the 8 scans in each on-off block. To do this, we set X_interest to a design matrix for the 8 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:

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

looks at the column of the design matrix for the 4th slice.

Analysing one run with fmrilm

Which_stats is a character matrix indicating which statistics for output, one row for each row of CONTRAST. If only one row is supplied, it is used for all contrasts. The statistics are indicated by strings, which can be anywhere in WHICH_STATS, and outputted in OUTPUT_FILE_BASEstring.ext, depending on the extension of INPUT_FILE. The strings are:

_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 all 16 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 effect (in %BOLD), try blob_brain:

Also of interest is the lag 1 autocorrelation smoothed by fwhm_cor=6.6168mm (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

The F statistics for testing for drift has df.F = 3 44. Note that we can use stat_threshold to find the threshold for significant drift (see later).

Thresholding the tstat image with stat_threshold and fdr_threshold

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, and spatial extent_threshold_1 is the extent of a single cluster chosen in advance.

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.

Producing an SPM style summary with stat_summary

Confidence regions for the spatial location of local maxima using conf_region

Extracting values from a file using extract

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

You also use blob_brain to colour code the t stat image thresholded at 4.93 with the delay in seconds:

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.

The more intersting question is how the connectivity might be modulated by the stimulus, that is, how the stimulus changes the connectivity. To do this, we must add an interaction (product) between the stimulus 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], i.e. 0 for the single covariate (reference voxel time course), then 1 to pick out the product of the stimulus with the covariate.

Unfortunately there is no evidence that the connectivity is modulated by the stimulus 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).