AFNI example

You will need the AFNI Matlab toolbox written by Ziad Saad.

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

Data contains all the fMRI data for the example if you want to duplicate the analysis.

The graphics are crude (orientation needs fixing) but of course you would use AFNI for that.

  • Looking at the fMRI data using pca_image
  • Making the design matrices using fmridesign
  • Analysing one run with fmrilm
  • Visualizing the results using view_slices, glass_brain and blob_brain
  • F-tests
  • Thresholding the tstat image with stat_threshold and fdr_threshold
  • Finding the exact resels of a search region with mask_resels
  • 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
  • Extracting values from a file using extract
  • Estimating the time course of the response
  • Estimating the delay
  • Effective connectivity of all voxels with a reference
  • Higher order autoregressive models

    Looking at the fMRI data using pca_image

    A Principal Components Analysis (PCA) of the fMRI data is sometimes useful as an exploratory tool that does not require any model for the data. It can reveal unexpected signals, or outlying frames. The PCA writes the data as a sum of temporal components times spatial components, chosen to best approximate the fMRI data. To define the brain, we use the first slice of the fMRI data as a mask. fmri_mask_thresh tries to find a nice threshold: input_file='c:/keith/kwafnidata/r1_time+orig.BRIK'; mask_file='c:/keith/kwafnidata/r1_time+orig.BRIK'; mask_thresh=fmri_mask_thresh(mask_file); pca_image(input_file, [], 4, mask_file, mask_thresh); [Click to enlarge 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):

    exclude=[1 2 3 4]; pca_image(input_file, exclude, 4, mask_file); [Click to enlarge image]

    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:

    X_remove=[ones(68,1) (1:68)']; X_interest=[zeros(4,8); repmat(eye(8),8,1)]; pca_image(input_file, exclude, 4, mask_file, [], [], [], X_remove, X_interest); [Click to enlarge image]

    Making the design matrices using fmridesign

    There are 68 scans, separated by TR=5 seconds, and 16 interleaved slicetimes which can be found from BrikInfo. The design was a block design with 4 scans off, 4 scans on, repeated 8 times, finishing with 4 scans off. The events are specified by a matrix of EVENTS, or directly by a stimulus design matrix S, or both. EVENTS is a matrix whose rows are events and whose columns are:

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

    [err, Info]=BrikInfo('c:/keith/kwafnidata/r1_time+orig.HEAD'); Info TR=5; frametimes=(0:67)*TR; slicetimes=Info.TAXIS_OFFSETS/1000; eventid=repmat(1,8,1); eventimes=((0:7)'*8+4)*TR; duration=ones(8,1)*4*TR; height=ones(8,1); events=[eventid eventimes duration height] X_cache=fmridesign(frametimes,slicetimes,events); plot(squeeze(X_cache.X(:,:,1,4)),'LineWidth',2) xlabel('frame number') ylabel('response') [Click to enlarge image]

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

    Analysing one run with fmrilm

    First set up the contrasts. CONTRAST is a matrix whose rows are contrasts for the statistic images, with row length equal to the number regressor variables in X_CACHE. EXCLUDE is a list of frames that should be excluded from the analysis. This must be used with Siemens EPI scans to remove the first few frames, which do not represent steady-state images.

    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.

    contrast=[1]; exclude=[1 2 3 4]; which_stats='_mag_t _mag_ef _mag_sd _cor _fwhm'; input_file='c:/keith/kwafnidata/r1_time+orig.BRIK'; output_file_base='c:/keith/results_test_afni/r1_time'; fmrilm(input_file, output_file_base, X_cache, contrast, exclude, which_stats); and the output gives fwhm_data = 6.6393 fwhm_cor = 6.6168 df = resid: 58 cor: 299 t: 52 The estimated fwhm of the data is fwhm_data=6.6393mm, and the amount of smoothing chosen for the autocorrelation parameters is fwhm_cor=6.6168mm. This achieves 52 df for the T statistic. The files generated are:

    c:/keith/results_test_afni/r1_time_cor+orig.BRIK c:/keith/results_test_afni/r1_time_fwhm+orig.BRIK c:/keith/results_test_afni/r1_time_mag_ef+orig.BRIK c:/keith/results_test_afni/r1_time_mag_sd+orig.BRIK c:/keith/results_test_afni/r1_time_mag_t+orig.BRIK plus header files.

    Visualizing the results using view_slices, glass_brain and blob_brain

    To look at slice 7 (zero based) or slice 8 (one based) of the t-statistic try

    t_file='c:/keith/results_test_afni/r1_time_mag_t+orig.BRIK'; m=fmris_read_image(t_file,8,1); imagesc(',[-6 6]); colorbar; axis xy; colormap(spectral); [Click to enlarge image]

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

    mask_file=input_file; clf; view_slices(t_file,mask_file,[],3,1,[-6 6]); [Click to enlarge image]

    or to see all 16 slices:

    clf; view_slices(t_file,mask_file,[],0:15,1,[-6 6]); [Click to enlarge image]

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

    glass_brain(t_file,3,mask_file); [Click to enlarge image]

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

    clf; blob_brain(t_file,5,mask_file); title('T>5, coloured by effect (%BOLD)'); [Click to enlarge image]

    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

    cor_file='c:/keith/test/results_afni/r1_time_cor+orig.BRIK'; clf; view_slices(cor_file,mask_file,0,7,1,[-0.15 0.35]); [Click to enlarge image]

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

    fwhm_file='c:/keith/test/results_afni/r1_time_fwhm+orig.BRIK'; clf view_slices(fwhm_file,mask_file,0,7,1,[0 20]); [Click to enlarge image]


    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 contrast.T=[0 1 0 0; 0 0 1 0; 0 0 0 1] which_stats='_mag_F'; output_file_base='c:/keith/test/results_afni/r1_time_drift'; fmrilm(input_file,output_file_base,X_cache,contrast,exclude,which_stats,cor_file); fwhm_data = 6.6393 fwhm_cor = 6.6168 df = resid: 58 cor: 299 t: [43 44 45] F: [3 44] clf; view_slices('c:/keith/test/results_afni/r1_time_drift_mag_F+orig.BRIK',mask_file,[],0:15,1,[0 50]); [Click to enlarge image]

    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

    Often we we want to define a search region by a mask. Mask_vol finds the volume of the search region, and the number of voxels of a mask_file above a mask_thresh - if no mask_thresh is supplied, it will attempt to fund one using fmri_mask_thresh. Using the raw fMRI data in mask_file to define the mask: [search_volume, num_voxels]=mask_vol(mask_file) stat_threshold(search_volume,num_voxels,6,52) peak_threshold = 4.9692 Cluster_threshold = 3.2549 peak_threshold_1 = 4.5133 extent_threshold = 412.5694 extent_threshold_1 = 118.1396 fdr_threshold(t_file,[],mask_file,[],52) fdr_thresh = 2.6985 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)

    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

    To look at all peaks and clusters above a threshold corresponding to an uncorrected P-value of 0.001 (default), together with a complete listing of P-values and false discovery rates (based on the local fwhm), use:
    stat_summary(t_file, fwhm_file, [], mask_file); [Click to enlarge image]

    clus vol resel p-val (one) 1 7580 22.16 0 ( 0) 2 5119 16.11 0 ( 0) 3 4725 12.65 0 ( 0) 4 1673 5.51 0 ( 0) 5 1673 5.18 0 ( 0) . . . 25 591 1.74 0.038 (0.001) 26 492 1.64 0.047 (0.001) 27 394 1.61 0.05 (0.001) 28 394 1.59 0.053 (0.001) 29 394 1.59 0.053 (0.001) . . . clus peak p-val (one) q-val (i j k) ( x y z ) 1 10.41 0 ( 0) 0 (32 26 7) ( -1.9 -20.6 3) 1 9.92 0 ( 0) 0 (30 28 7) ( 5.6 -13.1 3) 2 9.74 0 ( 0) 0 (30 30 1) ( 5.6 -5.6 45) 5 9.4 0 ( 0) 0 (37 30 0) (-20.6 -5.6 52) 5 9.09 0 ( 0) 0 (34 28 0) ( -9.4 -13.1 52) 2 8.76 0 ( 0) 0 (32 31 3) ( -1.9 -1.9 31) . . . 268 5.01 0.044 (0.012) 0 (33 43 5) ( -5.6 43.1 17) 115 4.99 0.046 (0.013) 0 (48 46 6) (-61.9 54.4 10) 257 4.97 0.05 (0.014) 0 (31 31 15) ( 1.9 -1.9 -53) 24 4.96 0.051 (0.014) 0 (40 30 4) (-31.9 -5.6 24) 20 4.91 0.061 (0.016) 0 (38 29 0) (-24.4 -9.4 52) . . .

    Confidence regions for the spatial location of local maxima using conf_region

    The spatial location of a peak also has a random error due to the noise component of the images. This error is roughly equal to FWHM / ( sqrt(4log(2)) * PEAK ), where FWHM is the effective fwhm of the data, and PEAK is the height of the peak. Note that higher peaks are more accurately located, as you would expect. This standard error assumes that the signal has a Gaussian spatial shape with a fwhm that matches FWHM. If this is not the case, a more robust method is based on the spatial shape of the peak itself. Instead of estimating the standard error, we can find a confidence region for the location by thresholding the peak at a value of sqrt( peak^2 - chisq_pvalue), where the probability of a chisq random variable with 3 df exceeding chisq_pvalue is the desired P-value (chisq_pvalue = 7.81 for P-value = 0.05). This is implemented by conf_region, which produces a 95% (P=0.05) confidence region image with values equal to the peak value in each confidence region. Thresholding this image at a value of say 4.91 in blob_brain, produces an image of confidence regions for locations of all peaks above 4.97. blob_brain colours the confidence regions by the peak height:

    conf_region(t_file,4.97,mask_file) conf_file='c:/keith/results_test_afni/r1_time_mag_t_95cr+orig.BRIK'; clf; view_slices(conf_file,mask_file,[],0:15) clf; blob_brain(conf_file,4.97,conf_file,4.97) title('Approx 95% confidence regions for peak location, coloured by peak height') [Click to enlarge image]

    Extracting values from a file using extract

    You may want to look at the effects and their standard deviations at some of the local maxima just found, for example the peak with value 10.46 in slice 7 (zero based) with voxel coordinates [32 26 7]:
    voxel=[32 26 7] ef=extract(voxel,'c:/keith/results_test_afni/r1_time_mag_ef+orig.BRIK') sd=extract(voxel,'c:/keith/results_test_afni/r1_time_mag_sd+orig.BRIK') ef/sd [df,spatial_av]=fmrilm(input_file,[],[],[],exclude); ref_data=squeeze(extract(voxel,input_file))./spatial_av*100; fitted=mean(ref_data)+ef*X_cache.X(:,1,1,voxel(3)+1); clf; plot(frametimes,[ref_data fitted],'LineWidth',2); xlim([15 max(frametimes)]); legend('Reference data','Fitted values'); xlabel('time (seconds)'); ylabel('fMRI response, percent'); title(['Observed (reference) and fitted data, ignoring trends, at voxel ' num2str(voxel)]); [Click to enlarge image]

    Estimating the time course of the response

    We may be intersted in estimating the time course of the response, rather than just the data as above. The easiest way to do this might be to average the responses beginning at the start of each `epoch' of the response, e.g. averaging the 8 scans starting at scans 1, 9, 17, ... . This does not remove drift, nor does it give you a standard error. A better, though more complex way, is to replace the block design of 1 block with say 8 blocks of 5s each covering the entire epoch, then estimate their effects with NO convolution by the hrf.

    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:

    eventid=kron(ones(8,1),(1:8)'); eventimes=((1:64)'+3)*TR; duration=ones(64,1)*TR; height=ones(64,1); events=[eventid eventimes duration height] X_bases=fmridesign(frametimes,slicetimes,events,[],zeros(1,5)); contrast=[eye(8)-ones(8)/8]; num2str(round(contrast*100)/100) which_stats='_mag_ef _mag_sd _mag_F'; output_file_base=['c:/keith/results_test_afni/r1_time01'; 'c:/keith/results_test_afni/r1_time02'; 'c:/keith/results_test_afni/r1_time03'; 'c:/keith/results_test_afni/r1_time04'; 'c:/keith/results_test_afni/r1_time05'; 'c:/keith/results_test_afni/r1_time06'; 'c:/keith/results_test_afni/r1_time07'; 'c:/keith/results_test_afni/r1_time08']; df=fmrilm(input_file,output_file_base,X_bases,contrast,exclude,which_stats,fwhm_rho) We need 12 output files, but only the effect and its standard error. The F-statistic has been requested so we can find voxels with significant responses. Its df's are df.F = 7 52, so its critical threshold is 7.38:

    stat_threshold(search_volume,num_voxels,6,[p df]) lm=locmax('c:/keith/results_test_afni/r1_time01_mag_F+orig.BRIK',7.38); num2str(lm) The extracted effects and their standard errors at the previous voxel are:

    values=extract(voxel,output_file_base,'_mag_ef+orig.BRIK') sd=extract(voxel,output_file_base,'_mag_sd+orig.BRIK') Plotting these out against time gives the estimated response, together with 1 sd error bars. Note that the times are offset by the slicetimes, so they are in seconds starting from the beginning of slice acquisition. On top of this we can add the modeled response by extracting the effects of stimulus, then multiplying these effects by the hrf convolved with the 20s block (also shown):

    time=(1:(8*TR*10))/10; X_hrf=fmridesign(time,0,[1 0 20 1]); hrf=squeeze(X_hrf.X(:,1,1)); plot((0:8)*TR+slicetimes(voxel(3)+1),values([1:8 1]),'k', ... [0:7; 0:7]*TR+slicetimes(voxel(3)+1), [values+sd; values-sd],'g', ... time,[ones(1,200) zeros(1,200)],'r', ... time,hrf*ef,'g','LineWidth',2); legend('Estimated response','Modeled response'); xlabel('time (seconds) from start of epoch'); ylabel('fMRI response, percent'); title(['Estimated and modeled response at voxel ' num2str(voxel)]); [Click to enlarge image]

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

    Estimating the delay

    Rather than the whole response at one voxel, we may prefer just one parameter, the delay, estimated at all voxels. Delays are shifts of the time origin of the HRF, measured in seconds. 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. Note that F statistics are not yet available with this option. If the length of CONTRAST_IS_DELAY is less then the number of contrasts, it is padded with zeros. Default is 0, i.e. all contrasts refer to magnitudes. Here is an example to estimate the delays of the stimulus for the above data. In this case every contrast is a delay, so CONTRAST_IS_DELAY = [1]:

    contrast=[1] which_stats='_del_t _del_ef _del_sd' output_file_base=['c:/keith/test/results_afni/r1_time'] fmrilm(input_file, output_file_base, X_cache, contrast, exclude, which_stats, cor_file); slice=7; subplot(2,2,1); view_slices('c:/keith/results_test_afni/r1_time_mag_t+orig.BRIK',mask_file,[],slice,1,[-6 6]); subplot(2,2,2); view_slices('c:/keith/results_test_afni/r1_time_del_ef+orig.BRIK',mask_file,[],slice,1,[-3 3]); subplot(2,2,3); view_slices('c:/keith/results_test_afni/r1_time_del_sd+orig.BRIK',mask_file,[],slice,1,[0 6]); subplot(2,2,4); view_slices('c:/keith/results_test_afni/r1_time_del_t+orig.BRIK',mask_file,[],slice,1,[-6 6]); [Click to enlarge image]

    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:

    clf; blob_brain('c:/keith/results_test_afni/r1_time_mag_t+orig.BRIK',5, ... 'c:/keith/results_test_afni/r1_time_del_ef+orig.BRIK'); title('Delay (secs) of hot stimulus where T > 5') [Click to enlarge image]

    Effective connectivity of all voxels with a reference voxel

    We might be interested in all voxels that are correlated with the time course at a pre-chosen reference voxel, removing the effect of the signal, to find those voxels that are effectively connected. To do this, we add the data from the pre-chosen voxel (ref_data) as a confound, then estimate its effect, sd, and t statistic as before. Note that ref_data is measured at times that depend on the slicetimes for that slice, i.e. slice 7 (zero based). For other slices, ref_data must be resampled at the same frametimes and slicetimes as the fmri data, so use fmri_interp to do this:
    output_file_base='c:/keith/test/results_afni/r1_time_connect'; contrast.C=1; which_stats='_mag_t _mag_ef _mag_sd' ref_times=frametimes'+slicetimes(voxel(3)+1); confounds=fmri_interp(ref_times,ref_data,frametimes,slicetimes); fmrilm(input_file, output_file_base, X_cache, contrast, exclude, which_stats, cor_file, [], confounds); clf; view_slices('c:/keith/test/results_afni/r1_time_connect_mag_t+orig.BRIK',mask_file,[],0:15,1,[-6 6]); [Click to enlarge image]

    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.

    output_file_base='c:/keith/test/results_afni/r1_time_connect_stim'; contrast.C=[0 1]; X_confounds=XC(X_cache,confounds); fmrilm(input_file, output_file_base, X_cache, contrast, exclude, ... which_stats, cor_file, [], X_confounds); clf; view_slices('c:/keith/test/results_afni/r1_time_connect_stim_mag_t+orig.BRIK',mask_file,[],0:15,1,[-6 6]); [Click to enlarge image]

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

    Higher order autoregressive models

    The last parameter of FMRILM is NUMLAGS, the number of lags or order of the autoregressive model for the temporal correlation structure. Order 1 (the default) seems to be adequate for 1.5T data, but higher order models may be needed for 3T data. To assess this, the first NUMLAGS autoregression parameters can be outputted as 'frames' by setting the 8th element of WHICH_STATS to 1. These parameters should be zero for orders higher than the true order of the data. Note: if NUMLAGS>1 it takes a lot longer to fit. Here is an example of an AR(4) model (note that setting the intervening parameters to [] uses their defaults): which_stats='_mag_t _mag_ef _mag_sd _cor _AR' contrast=[1]; output_file_base='c:/keith/test/results_afni/r1_time_ar4'; fmrilm(input_file, output_file_base, X_cache, contrast, exclude, ... which_stats, [], [], [], [], [], [], 4); fwhm_data = 6.6393 fwhm_cor = 12.9939 df = resid: 58 cor: 1478 t: 52 Note that the amount of smoothing (fwhm_cor = 12.9939) is more than for the AR(1) model (fwhm_cor = 6.6168) because AR(4) models have more parameters, so more smoothing is needed to get the same df.t = 52.

    clf; view_slices('c:/keith/test/results_afni/r1_time_ar4_AR+orig.BRIK',mask_file,0,7,1:4,[-0.15 0.35]); [Click to enlarge image]

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