All matlab programs are zipped or in the
toolbox directory (previous version here)
- you need to copy ALL the files to be able to read and write AFNI, ANALYZE and MINC formatted files.

If you are using
MINC formatted files you will need the extra tools in the
emma toolbox also available for
Windows and Linux.
Installing the Windows and Linux versions is a bit tricky.
In the zipped file that you download, there are .dll files and .exe files -
these must be moved to a directory in the path, e.g. matlab/bin/win32 directory.
Also you must create a directory c:/tmp or /tmp for
emma to use as scratch space.

Summary of method

The fMRI data was first converted to percentage of whole volume.
The statistical analysis of the percentages was based on a linear model
with correlated errors. The design matrix of the linear model was first
convolved with a hemodynamic response function modelled as a difference of
two gamma functions timed to coincide with the acquisition of each slice.
Temporal drift was removed by adding a cubic spline in the frame times
to the design matrix (one covariate per 2 minutes of scan time),
and spatial drift was removed
by adding a covariate in the whole volume average. The correlation
structure was modelled 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,
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 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.

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. By default, the data is normalised at each voxel by
dividing by the sd over time, so that the PCA is done in the spatial correlation matrix
rather than the spatial covariance matrix (setting the parameter ISCOV=1 will
give the unnormalised, covariance matrix analysis). 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, in this case 452.
Pca_image produces the following
graph of the first 4 components:
input_file='c:/keith/data/test_19971124_1_100326_mri_MC.mnc';
mask_file=input_file;
mask_thresh=fmri_mask_thresh(mask_file);
pca_image(input_file, [], 4, mask_file, mask_thresh);
colormap(spectral);

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

exclude=[1 2 3];
V=pca_image(input_file, exclude, 4, mask_file);

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.

The hrf is modeled as the difference of two gamma density functions (Glover, G.H. (1999).
"Deconvolution of impulse response in event-related
BOLD fMRI." NeuroImage, 9:416-429).
The parameters of the hrf are specified by a row vector whose elements are:

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:

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:

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.
We wish to look for regions responding to the
hot stimulus, the warm stimulus, and the difference hot-warm:

contrast=[1 0;
0 1;
1 -1];
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. Default is [1], but [1 2 3] is better:

exclude=[1 2 3];
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.

which_stats='_mag_t _mag_sd _mag_ef _mag_F _cor _fwhm';
Input minc file for one session, and output base for each row of CONTRAST
(hot for hot, wrm for warm, hmw for hot-warm):

fmrilm(input_file, output_file_base, X_cache, contrast, exclude, which_stats)
and the output gives
fwhm_data =
8.7855
fwhm_cor =
12.4444
df =
resid: 110
cor: 1235
t: [103 103 99]
F: [2 102]
The estimated fwhm of the data is fwhm_data=8.7855mm, and the amount of smoothing chosen for the
autocorrelation parameters is fwhm_cor=12.4444mm. This achieves df's of 103 103 99 for each of the T statistics, and 2 102 df for the F statistic. Note that the numerator df of F is 2, not 3, because
the third contrast is a linear combination of the first two. The files generated are:

F-tests

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 temproal trend or drift, since we would be interested in
either a linear, quadratic or cubic trend, or any linear combination
of these. The drift terms appear in the order
constant, linear, quadratic, cubic etc. so to test for a drift, use 'contrast.T'
(T for temporal trend) as follows:
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/test_100326_drift';
fmrilm(input_file,output_file_base,X_cache,contrast,exclude,which_stats,cor_file)
fwhm_data =
8.7855
fwhm_cor =
12.4444
df =
resid: 110
cor: 1235
t: [95 95 96]
F: [3 95]
Note that by setting the fwhm_cor parameter with the cor_file name
c:/keith/test/results/test_100326_hot_cor.mnc we avoid re-calculating the autocorrelations.
To visualize the regions with significant drift:
clf;
view_slices('c:/keith/test/results/test_100326_drift_mag_F.mnc',mask_file,[],0:11,1,[0 50]);
stat_threshold(1000000,26000,8,[3 95])

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

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

events=[zeros(20,1)+1 eventimes duration ones(20,1);
zeros(20,1)+2 eventimes duration temperature;
zeros(20,1)+3 eventimes duration temperature.^2]
contrast=[0 0 1];
The contrast will test for a quadratic effect. Taking this further, we may wish
to test for a cubic (Figure c), quartic (Figure d) or even a 'non-linear' effect
of temperature, where the temperature effects are arbitrary. Note that the
quartic model is identical to the arbitrary model because a quartic can be
fitted exactly to any arbitrary 5 points. To model the arbitrary effect,
simply assign a different event type from 1 to 5 for each of the 5 different
temperature values:
events=[floor((temperature-35)/3.5)+1 eventimes duration ones(20,1)]
To look for any effect of the stimulus compared to baseline, use:
contrast=[1 0 0 0 0;
0 1 0 0 0;
0 0 1 0 0;
0 0 0 1 0;
0 0 0 0 1]
contrast=eye(5)
which_stats='mag_F';
The resulting F statistic image will have 5 and 107 d.f.. To see if the
changes in temperature have any effect, subtract the average from each
row of the contrasts:
contrast=[.8 -.2 -.2 -.2 -.2;
-.2 .8 -.2 -.2 -.2;
-.2 -.2 .8 -.2 -.2;
-.2 -.2 -.2 .8 -.2;
-.2 -.2 -.2 -.2 .8]
contrast=eye(5)-ones(5)/5
which_stats='mag_F';
The resulting F statistic image will have 4 and 107 d.f.. To see if the
changes in temperature have a linear effect, subtract the average
from the 5 temperature values:
contrast=[-7.0 -3.5 0 3.5 7.0];
which_stats='_mag_t _mag_ef _mag_sd';
This will produce a T statistic with 107 d.f.. How is this different
from the previous linear effect? The effect is the same, but the
standard deviation may be slightly different. The reason is that in
this model, we are basing the standard deviation on errors about each
temperature value; in the first analysis, it is based on errors about
a linear model. As a result, the degrees of freedom is slightly
different: 112 before, but 109 here.

multistat(input_files_ef,input_files_sd,[],[],X,contrast,output_file_base,which_stats,Inf)
Note that the program prints and returns the final
degrees of freedom of the tstat image, which is 396 (see discussion
later):

df =
resid: 3
fixed: 396
sdratio: Inf
t: 396
For a more elaborate analysis, e.g. comparing the first two sessions with
the next two, you can do it by:
X=[1 1 0 0; 0 0 1 1]';
contrast=[1 -1];

Fixed and random effects

In the above run, we did a fixed effect analysis, that is, the
analysis is only valid for the particular group of 4 sessions on the test subject. It is not valid for an arbitrary session on the test subject, nor for
the population of subjects in general. It may turn out that the
underlying effects themselves vary from session to session, and from
subject to subject. To take this extra source of variability into
account, we should do a random effects analysis.

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,

Suppose we search c:/keith/test/results/test_multi_hmw_t.mnc for local maxima
inside a search region of 1000cc.
The voxel volume is 2.34375*2.34375*7= 38.4521, so the number of voxels
is 1000000/38.4521=26000.
There is roughly 8 mm smoothing.
The degrees of freedom returned by multistat,
is 100. By default the significance is the standard 0.05. The threshold to use is

stat_threshold(1000000,26000,8,100)
which gives 4.89. 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:
mask_file='c:/keith/data/test_19971124_1_100326_mri_MC.mnc';
[search_volume, num_voxels]=mask_vol(mask_file)
gives search_volume = 1183800 mm^3 (just over a litre) and num_voxels = 30786.
Using these values,
stat_threshold(search_volume,num_voxels,8,100)
gives a threshold of 4.93.

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:

stat_threshold(search_volume,num_voxels,8,[100; 100])
The thresholds for peak height are the same, but the
thresholds for spatial extent increase slightly to extent_threshold =
683 mm^3 and extent_threshold_1 = 223 mm^3 for a pre-chosen cluster.
If the FWHM comes from
multistat, two degrees of freedom are required: the first is the
residual degrees of freedom df.resid, equal to the
number of runs - number of parameters in the model = 4 - 1 = 3, and the second is df.t, the effective
degrees of freedom after smoothing the rfx = 100 (output as the
first parameter of multistat). These are put into the second row of a 2 x 2 matrix as follows
(the 0 indicates that we want calculations for a t statistic):
stat_threshold(search_volume, num_voxels,8,[100 0; 3 100])
The thresholds for local maxima are the same, but the thresholds for spatial extent
jump to 1020 mm^3 and 239 mm^3 for a pre-chosen cluster.

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:

t_file='c:/keith/test/results/test_multi_hmw_t.mnc';
fdr_threshold(t_file,[],mask_file,[],100)
The result is a threshold of 2.66, considerably lower than the random field /
Bonferroni threshold of 4.93 above. But remember that the interpretation is
different: at least 95% of the voxels above 2.66 are true signal (on average), whereas the
19 times out of 20, there will be *no* false positives in the voxels above 4.93.

For the perfectionist, there is a way of finding the exact resels
of a search region in dimensions 0,1,2,3. These are used by
stat_threshold
to get more accurate thresholds and P-values that take into account
the shape of the search region, not just its volume (see Worsley et al., 1996, above).
Up to now, stat_threshold has been assuming that
the search region is a ball, so that the resulting thresholds are slighty too
low (liberal) for non-ball shaped regions.
Mask_resels will find the exact resels of
any search region defined by a mask_file above mask_thresh (if mask_thresh is not
provided, it will use fmri_mask_thresh
to find one). It works by filling the
mask with a tetrahedral mesh that is warped so it touches the boundary, thus
overcoming the 'jagged boundary' problem of the method in Worsley et al. (1996).
The FWHM must also be supplied, here assumed to be 8mm:
mask_file='c:/keith/data/test_19971124_1_100326_mri_MC.mnc';
resels=mask_resels(8,[],mask_file);
the result is resels = [1 36.3 516.1 2291.6]. The first component (dimension 0)
is the Euler characteristic of the mask, which is 1 if the search region is simply
connected (topologically equivalent to a ball). The second component (dimension 1)
is 2 * caliper diameter / FWHM, the third (dimension 2) is 1/2 * surface area / FWHM^2,
and the fourth (dimension 3) is the familiar volume / FWHM^3. Note that this is
not quite the same as search_volume / 8^3 = 1183800 / 8^3 = 2312.1 found using
mask_vol. Mask_resels is more accurate - it fills the mask with tetrahedra out to
the boundary (interpolated between voxels, a bit like 'marching cubes'), whereas
mask_vol assumes that the search region is a set of box-shaped voxels surrounding
each voxel in the search region (so it has a 'jagged' boundary). Putting these
resels into stat_threshold in place of search_volume, and setting FWHM=1:
stat_threshold(resels,num_voxels,1,100)
gives peak_threshold = 4.93, peak_threshold_1 = 4.30 and
extent_threshold = 1.2425, extent_threshold_1 = 0.4310
(note these are now measured in resels, not mm^3). In other words, the results
are almost identical to what we had before. The reason is that
stat_threshold
reports the lower of the Bonferroni
threshold and the random field threshold, and in this case the Bonferroni
threshold, which depends on the number of voxels and not the resels, is lower.
To see the
difference in the random field threshold, we can set the number of voxels to infinity:
stat_threshold(resels,Inf,1,100)
5.2162
stat_threshold(search_volume,Inf,8,100)
5.2193
The returned thresholds are almost identical, so again this is just for the
perfectionist!

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:

which_stats=which_stats='_wresid';
contrast=[];
output_file_base='c:/keith/test/results/test_100326_hot';
fmrilm(input_file, output_file_base, X_cache, contrast, exclude, which_stats);
fwhm_info='c:/keith/test/results/test_100326_hot_wresid.mnc';
resels=mask_resels(fwhm_info,[],mask_file);
The result is resels =[ 1 -476.6 304.0 1639.8] which has considerable negative
caliper diameter for reasons that are not clear at the moment. Using this
again in stat_threshold,
and allowing for the variability in estimating the resels:
stat_threshold(resels,num_voxels,1,[100; 100])
The peak thresolds are the same as before because again the Bonferroni wins, but
the extent_threshold = 1.2108 and extent_threshold_1 = 0.4349 are a bit smaller
due to the smaller 3D resels. Note that the
extents are in resels, not mm^3, which allows for non-uniform FWHM.
It would not be complete without an example from multistat:
which_stats='_wresid';
contrast=[];
output_file_base='c:/keith/test/results/test_multi_hmw';
multistat(input_files_ef,input_files_sd,[],[],X,contrast,output_file_base,which_stats);
fwhm_info='c:/keith/test/results/test_multi_hmw_wresid.mnc';
resels=mask_resels(fwhm_info,[],mask_file)
stat_threshold(resels,num_voxels,1,[100 0; 3 100])
The result is resels =[1 95.4 311.7 1561.4 ], similar to what we had before, with
identical peak thresholds (because Bonferroni wins) and
extent_threshold = 1.7204 and extent_threshold_1 = 0.4669 presumably
because there is more uncertainty in estimating the resels.

To locate peaks and clusters masked by the first slice if the fMRI
data, use
lm=locmax(t_file, 4.94, mask_file)
which gives a list of local maxima above 4.94, sorted in descending order.
The first column of LM is the values of the local maxima, and the next three
are the x,y,z voxel coordinates starting at 1.
The 5th column is a cluster id, and the 6th is the cluster volume.
To find P-values for these local maxima, add the list of
peak values from the first column of lm, e.g.
pval=stat_threshold(1000000,26000,8,100,lm(:,1))
Locmax can now produce a file
that has the cluster id as its values, useful for locating the clusters.
It can also produce cluster
resels as well as volumes if a FWHM scalar or file is given - see the help.

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

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 9 in register or
blob_brain,
produces an image of confidence regions for locations of all peaks above
9. blob_brain colours the confidence regions
by the peak height:

We might be interested in the common regions that are activated by both the
hot and the warm stimuli, called a conjunction. In general, we might be interested in the common regions activated by an arbitrary number of arbitrary test statistic images. These images may be the same contrast on different runs/subjects, or different contrasts from the same run. The test statistics may have different df, or even be a mixture of T and F statistics. The best way to do this is to first convert all test statistic images to P-values using stat_summary as above. Then the conjunction can be obtained by taking the maximum of the p-values (or minimum of the negative p-values) and thresholding at the desired overall p-value of say 0.05. The surviving voxels are guaranteed to show activation on all the test statistic images, even if they are correlated. The test can be quite stringent, since if just one image is unactivated we will fail to detect the conjunction.

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

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.6 in slice 3 (zero based) with voxel coordinates [62 69 3]:
voxel=[62 69 3]
ef=extract(voxel,'c:/keith/test/results/test_multi_hmw_ef.mnc')
sd=extract(voxel,'c:/keith/test/results/test_multi_hmw_sd.mnc')
ef/sd
ef =
1.1198
sd =
0.1047
ans =
10.6969
which gives 1.1198 and 0.1047 respectively. Note that 1.1198/0.1047 = 10.7,
the value of the T statistic. A matrix of voxels can be given, and the corresponding
vector of values will be extracted.

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

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 12 scans starting at scans 1, 13, 25, ... 109.
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 2 blocks (9s hot, 9s warm) with
say 12 blocks of 3s each covering the entire epoch, then estimate their
effects with NO convolution by the hrf:
eventid=kron(ones(10,1),(1:12)');
eventimes=frametimes';
duration=ones(120,1)*3;
height=ones(120,1);
events=[eventid eventimes duration height]
events =
1 0 3 1
2 3 3 1
3 6 3 1
4 9 3 1
5 12 3 1
6 15 3 1
7 18 3 1
8 21 3 1
9 24 3 1
10 27 3 1
11 30 3 1
12 33 3 1
1 36 3 1
2 39 3 1
3 42 3 1
.
.
.
X_bases=fmridesign(frametimes,slicetimes,events,[],zeros(1,5));
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:

stat_threshold(search_volume,num_voxels,8,df.F)
which gives 5.27. The local maxima above this are:

lm=locmax([output_file_base(1,:) '_mag_F.mnc'],5.27);
num2str(lm)
There is a peak of 9.24 at [61 67 3], very near the previous voxel [62 69 3].
The extracted effects and their standard errors at the previous
voxel are:

voxel=[62 69 3]
values=extract(voxel,output_file_base,'_mag_ef.mnc')
sd=extract(voxel,output_file_base,'_mag_sd.mnc')
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 the hot and warm stimuli, then
multiplying these effects by the hrf convolved with the 9s
blocks for hot and warm stimuli (also shown):

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

which_stats='_del_ef _del_sd _del_t _mag_t';
then the magnitude images are not quite the same as what you would get from fmrilm
without asking for delays, i.e.
which_stats='_mag_t';
The reason is that in order to calculate delays, fmrilm
replaces the convolution of the stimulus with the hrf by
the convolution of the stimulus with
two basis functions for the hrf, in the linear model. It uses
these two basis functions for every response that has a non-zero
value in any contrast in the delays, otherwise it just uses one basis
function (the hrf itself). If you want to force two basis functions
into the model for any extra responses, even though you don't want to estimate
delays in these responses, see the next parameter
NUM_HRF_BASES.

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

register c:/keith/test/results/test_100326_hot_mag_t.mnc c:/keith/test/results/test_100326_hot_del_ef.mnc
After synchronising the two images,
you can then search the tstat image in the first window
for regions of high signal e.g. >5, then
read off the estimated delay at that location from the second window.
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/test/results/test_100326_hot_mag_t.mnc',5, ...
'c:/keith/test/results/test_100326_hot_del_ef.mnc');
title('Delay (secs) of hot stimulus where T > 5')

The efficiency of a design for an effect is inversely proportional to the sd of that effect
using such a design.
The program efficiency calculates the sd of an effect
without needing any data. It simply executes all the code of fmrilm that
calculates the sd, setting the sd of the noise to 1, so the returned sd
is relative to the sd of the noise. The call
is a lot like fmrilm except that you don't supply input_file or
output_file_base, but you must supply details of the design through
X_cache from fmridesign, the contrast, and the frames to exclude.
Because it works out the sd at a single voxel, it is quite fast. It returns a matrix of sd's
for each contrast (rows) and slice (columns). For the design above:
efficiency(X_cache, contrast, exclude)
ans =
Columns 1 through 8
0.1558 0.1565 0.1559 0.1566 0.1560 0.1567 0.1561 0.1567
0.1619 0.1618 0.1619 0.1617 0.1619 0.1617 0.1618 0.1616
0.1918 0.1916 0.1918 0.1916 0.1918 0.1916 0.1917 0.1915
Columns 9 through 13
0.1562 0.1567 0.1563 0.1567 0.1564
0.1618 0.1615 0.1618 0.1613 0.1618
0.1917 0.1915 0.1917 0.1914 0.1917
The sd of the hot stimulus is about 0.156 in all 13 slices,
the warm stimulus is about the same (0.162) and the
sd for the hot-warm stimulus is larger (0.192) by about sqrt(2) since it is the difference
of two effects.

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

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 functionally 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 3 (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.
Note also how the "C" at the end of the "contrast" referes to the confounds, so a value
of 1 picks out the confound we just added:
output_file_base='c:/keith/test/results/test_100326_connect';
contrast.C=1;
which_stats='_mag_ef _mag_sd _mag_t';
voxel=[62 69 3];
ref_times=frametimes'+slicetimes(voxel(3)+1);
[df,spatial_av]=fmrilm(input_file,[],[],[],exclude);
ref_data=squeeze(extract(voxel,input_file))./spatial_av*100;
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/test_100326_connect_mag_t.mnc',mask_file,[],0:11,1,[-6 6]);

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 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_ef _mag_sd _mag_t _cor _AR';
contrast=[1 -1];
output_file_base='c:/keith/test/results/test_100326_hmw_ar4';
fmrilm(input_file, output_file_base, X_cache, contrast, exclude, ...
which_stats, [], [], [], [], [], [], 4);
fwhm_data =
8.7855
fwhm_cor =
14.5716
df =
resid: 110
cor: 1824
t: 99
Note that the amount of smoothing (fwhm_cor = 14.5716) is
more than for the AR(1) model (fwhm_cor = 12.4444) because AR(4) models have more
parameters, so more smoothing is needed to get the same df.t = 99.

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