**Random Effects Analyses with Multiple Basis
Functions per subject**

This page is available as a WORD file: spm2_face_rfx.doc

*Rik Henson mailto:r.henson@fil.ion.ucl.ac.uk*

*Wellcome** Department*

*Institute
of Cognitive Neuroscience*

*Will
Penny mailto:w.penny@fil.ion.ucl.ac.uk*

*Wellcome** Department*

*May
2004*

CONTENTS

Example 1.1: Main Effect: Canonical HRF

Example 1.2: Differential Effect: Canonical HRF

Example 2.1: Main Effect: Informed
basis set

Effects of Interest F-contrast

Example 2.2: Differential Effect: Informed basis set

Example 3.1: Main Effect: FIR Model

Example 3.2: Differential Effect: FIR Model

These examples illustrate multisubject "random effects" in SPM2.

They illustrate 2nd-level ("random effects") analyses of event-related responses, as characterised by one or more temporal basis functions, across 12 subjects.

The examples consist of three basic types of 2nd-level model (M2):

**M2c**. Using contrast images for the
canonical HRF only

1 observation (contrast image) per subject

(one-sample t-test)

**M2i**. Using contrast images for the
"informed" basis set, consisting of the canonical HRF and its two
partial derivatives with respect to time (onset latency) and dispersion

3 observations (contrast images) per subject

(one-way 1x3 ANOVA)

**M2f**. Using contrast images from a very
general "Finite Impulse Response" (FIR) basis set, with 12 x 2 second
timebins

12 observations (contrast images) per subject

(one-way 1x12 ANOVA)

M2c is the same model used to illustrate parametric "random effects" in SPM99 and nonparametric "random effects" in SNPM in the example dataset: http:/www.fil.ion.ucl.ac.uk/spm/data/#SPM00Multi (though the results differ slightly because the data - ie contrast images - derive from slightly different first-level models).

M2i and M2f are used to illustrate how SPM2 deals with nonsphericity.

M2i can also be used to illustrate the sufficiency of SPM's "informed" basis set (ie canonical HRF and two derivatives) to capture the majority of event-related variance in the brain (at least within these data).

Two types of 1st-level contrasts are tested with each type of 2nd-level model (see below), resulting in 6 example SPM analyses.

The data come from the `implicit' condition of the study:

Henson, R.N.A, Shallice, T., Gorno-Tempini,
M.-L. & Dolan, R.J (2002). Face repetition effects in implicit and explicit
memory tests as measured by fMRI. *Cerebral
Cortex*, 12, 178-186

although the 1st-level design matrices (and therefore resulting contrast images) used do not correspond exactly to those used in that study.

Unlike the single-subject Fixed Effects example dataset, only two event-types were modelled: famous and nonfamous faces (initial and repeated presentations were collapsed together, as were correct and incorrect responses). Briefly, greyscale photographs of 52 famous and 52 nonfamous face were presented for 0.5s for fame judgment task (one of two right finger key presses). The minimal SOA (SOAmin) was 4.5s, with all faces randomly intermixed together with a further 52 null events (ie 2/3 probability of a face every SOAmin).

Original images were continuous EPI (TE=40ms,TR=2s) 24 descending slices (64x64 3x3mm2), 3mm thick, 1.5mm gap.

2nd-level models M2c and M2i derive from a 1st-level model (M1i), in which the events were modelled with Nf=3 basis functions: the canonical HRF, its partial derivative with respect to onset latency ("temporal derivative") and its partial derivative with respect

to dispersion ("dispersion derivative").

2nd-level model M2f derives from an alternative 1st-level model (M1f), in which the same events were modelled with Nf=12 basis functions instead: corresponding to 2s timebins from 0-24s poststimulus (SPM's "Finite Impulse Response" or FIR basis set).

In both first-level models (M1i and M1f), the contrast images (con*.img's) come from session-specific contrasts within a large (multisession) 1st-level Fixed Effects design matrix, with one session per subject. (Note that the resulting con*.img's could equally well have been produced from 12 separate 1st-level models, one per subject.)

For each type of model, two types of 1st-level contrast are examined:

1. The main effect of faces versus baseline (eg, a [0.5 ... 0.5] contrast for each basis function, or "kron(eye(Nf),[0.5 0.5])" more generally)

2. The differential effect of famous versus nonfamous faces (eg, a [-1 ... 1] contrast for each basis function, or "kron(eye(Nf),[-1 1])" more generally)

The 12 (subjects) x 3 (basis functions) x 2 (contrast-types) con*.imgs from the 1st-level model using the informed basis (M1i) set are in the zipped file cons_informed.zip

The 12 (subjects) x 12 (basis functions) x 2 (contrast-types) con*.imgs from the 1st-level model using the FIR basis (M1f) set are in the zipped file cons_fir.zip

Each contrast-type is examined in a separate SPM analysis (resulting in two SPM analyses per type of second-level model).

In this example, only contrasts that involve the canonical HRF basis function in the 1st-level model M1i are taken to a second-level analysis.

For the main effect versus baseline, these happen to correspond to the contrast images numbered 3-14 in 1st-level model M1i, ie:

con_0003.img (canonical HRF, subject 1)

con_0004.img (canonical HRF, subject 2)

...

con_0014.img (canonical HRF, subject 12)

(available here: cons_informed.zip)

These images comprise the data for M2c, which is simply a "one-sample t-test" in SPM2.

Type SPM at the matlab prompt.

Now change to a new directory (this is easy
to forget !)

Press the 'Basic models' button.

Select design type ..... [One sample
t-test]

Then select the 12 images (con_0003.img to
con_0014.img)

GMsca: grand mean scaling [No]

explicitly mask images [No]

Global calculation [omit]

SPM will then show you the design matrix (simply a single column of 1's which will appear as a white box on a white background). This design is encoded in the "SPM.mat" file that is written to the currect directory.

Then press the 'Estimate' button, and
select the "SPM.mat" file.

SPM will now estimate the parameters (ie. the size of the population effect at each voxel - simply the average of the con*.img's in this model).

Now press the 'Results' button.

Select the SPM.mat file.

In the contrast manager press 'Define new
contrast' (select T).

Enter [1] in the contrast section and enter
'activation' as a 'name'.

A [1] contrast tests for "activations" vs baseline, a [-1] for "deactivations".

Press the '..submit' button.

Press OK.

Now press the 'Done' button.

Mask with other contrast(s) [No]

Title for comparison [activation: Can HRF
only]

p value adjustment to control [FWE]

Family-wise p-value [0.05]

& Extent threshold

SPM will now display the thresholded t-statistic image. This shows the voxels that are significantly active (correcting for multiple comparisons across all voxels) in the population from which the subjects were drawn. They include bilateral posterior fusiform (e.g, 30 -63 -27, Z=6.15), SMA, and, at a more liberal threshold, left motor cortex).

The contrast images for the difference between famous and nonfamous faces happen to correspond to contrast numbers 39-50 in the 1st-level model M1i, ie:

con_0039.img (canonical HRF, subject 1)

con_0040.img (canonical HRF, subject 2)

...

con_0050.img (canonical HRF, subject 12)

The model ("one-sample t-test") is the same as in the main effect analysis in Eg1.1 above (so repeat the same steps, except selecting con_0039.img to con_0050.img, instead of con_0003.img to con_0014.img).

If you then evaluate the T-contrasts [1] and [-1], you will reveal voxels in which "Famous > Nonfamous" and "Nonfamous > Famous" respectively (at least, in terms of the loading the canonical HRF), since the contrasts at the first-level that created these contrast images were "Famous - Nonfamous". (You would get equivalent results if your contrast images were "Nonfamous - Famous" instead, and you reversed the contrast weights at the 2nd-level.) The F-contrast [1] is the two-tailed equivalent (Famous <> Nonfamous).

Note that there are only a few voxels (small clusters) surviving a FWE corrected p value of 0.05 in the [1] contrast (but cf. below); and no voxels surviving correction in the [-1] contrast.

For this example, 3 contrast images per subject are taken to the 2nd-level.

For the main effect versus baseline, these happen to correspond to the contrast images numbered 3-14 in the 1st-level model M1i, ie:

con_0003.img (canonical HRF, subject 1)

con_0004.img (canonical HRF, subject 2)

...

con_0015.img (temporal derivative, subject 1)

con_0016.img (temporal derivative, subject 2)

...

con_0027.img (dispersion derivative, subject 1)

con_0028.img (dispersion derivative, subject 2)

...

These images comprise the data for M2i, which is a "one-way ANOVA" in SPM2 (no constant term).

Note that no constant term should be included in the model. If one added a constant term, the contrasts in the 2nd-level model would have to sum to zero. This would entail contrasts that compared one basis function with another (e.g, F-contrasts [1 -1 0; 0 1 -1; -1 0 1], which test the null hypothesis that the mean of the basis functions differ). This is not what we want (and is pretty meaningless). We want to ask if the means of one or more of the basis functions are different from zero (e.g, F-contrasts [1 0 0; 0 1 0; 0 0 1] - see below - which test the null hypothesis that the three means are zero).

The three basis functions (in the 1st-level model) are scaled differently (more precisely, their relative scaling is arbitrary). This means that their parameter estimates are likely to have different variances over subjects. This is called inhomogeniety of variance, and invalidates classical parametric statistics, unless some correction is made. It is one example of nonsphericity, which SPM2 "accommodates" by estimating the error covariance (using ReML) and using it to "pre-whiten" the data and design matrix (see Glaser et al. (2003) Variance Components, Chapter 39, Human Brain Function, 2nd Edition, Elsevier for further details, available here Ch9.pdf). This facility is a major extension over SPM99, and is better for multi-subject designs with more than one parameter per subject (sphericity was assumed in SPM99, which is unlikely to be true

for repeated measures designs).

The values of the parameter estimates for each basis function are likely to be correlated across subjects (since they are "repeated measures"). The fact that the functions themselves are orthogonal is irrelevant here: the correlation in the parameter estimates is simply because, if one subject has a larger loading (parameter estimate) on the canonical HRF, they are also likely to have larger loadings on the derivatives (if their response differs from the canonical). These correlations induce nonzero covariances in the error.

This is another example of nonsphericity, which SPM2 "accommodates" with extra hyper-parameters in the manner described above (see also Glaser et al, 2003, Ch9.pdf).

Type SPM at the matlab prompt.

Now change to a new directory (this is easy
to forget !)

Press the 'Basic models' button.

Select design type ..... [One way ANOVA]

#group's ... enter [3]

group 1: select images: select con_0003.img
to con_0014.img

group 2: select images: select con_0015.img
to con_0026.img

group 3: select images: select con_0027.img
to con_0038.img

GMsca: grand mean scaling [No]

Threshold masking [none]

explicitly mask images [No]

Global calculation [omit]

non-sphericity correction [yes]

replications are over?... [repl(12)]

correlated repeated measures [yes]

The "yes" to the first "non-sphericity?" question allows for inhomogeniety of variance (ie 3 separate hyper-parameters over the "block diagonal" terms of the error covariance matrix).

The "yes" to the "repeated measures?" question allows for inhomogeniety of covariances (ie 3 further hyper-parameters over the three "off-diagonal" block terms of the error covariance matrix).

The "repl(12)" answer to the "replications are over?" question tells SPM2 that the replications come from the (implicit) factor of "replications", which in this case we know corresponds to subjects (and the other factor, group (=basis function), is the repeated measure).

The answers to these questions result in
the specification of a set of matrices (bases) that characterise the covariance
matrix. The first three correspond to the variance of each of the canonical,
temporal and dispersion derivatives:
SPM.xVi.Vi

The next three correspond to covariances: SPM.xVi.Vi

derivatives).

After estimation the actual covariance values (hyper-parameters) are given by SPM.xVi.h (the six entries correspond to the above bases). Note that these are 'global' values which are scaled by a voxel specific-value to achieve a model covariance that best matches the empirical covariance at each voxel.

After estimation, test the "effects of interest" F-contrast:

In the contrast manager press 'Define
new contrast' (select F).

Enter ["eye(3)"] (which in matlab
evaluates to [1 0 0; 0 1 0; 0 0 1])

as the contrast weights, and
"Can+Tem+Dis" as the name.

This contrast will reveal voxels that show some form of event-related response that can be captured by (ie, lies in the space spanned by) the three basis functions (e.g, 30 -60 -27, Z=7.43).

Note also how the design matrix changes after estimation. This is because it has been pre-whitened. In particular, the (barely visible) off-diagonal entries in the design matrix give an indication of the degree of correlation between the basis functions across subjects.

Note that because the data have also been pre-whitened our interpretation of the parameter estimates (the 'betas') is unchanged. Effectively the parameters have been estimated using 'Weighted Least Squares (WLS)', where the weights relate to the estimated error covariance structure. SPM2 implements WLS by pre-whitening the data and the design matrix and then using 'Ordinary Least Squares' (OLS).

It is also informative to evaluate the
T-contrast [1 0 0]

(ie positive loadings on the canonical HRF
only).

At a FWE correct p-value of 0.05, note more voxels (including now left motor cortex) and higher Z-values (e.g, 39 -57 -30, Z=7.53) for this main effect vs baseline compared to the equivalent contrast ([1]) in the model that uses only the canonical HRF (in Eg1.1 above). The main reason for this increased power is the increase in the degrees of freedom, which entails better estimators of the underlying error (co)variance. The price of this increased power is a stronger assumption about the nonsphericity, namely that it has the same structure across (activated) voxels - the "pooling device" (see Glaser et al. (2003) Ch9.pdf).

Finally, evaluate the F-contrasts [0 1
0] and [0 0 1].

These contrasts reveal voxels that load (positively or negatively) on the temporal and dispersion derivatives respectively. These contrasts reveal that there is significant variability (at p<.05 corrected) that is not captured by the canonical HRF alone (see Eg3.1 below for more discussion; see also to Henson et al (2000): hbm-fir.pdf).

In other words, some regions have earlier or later, or wider or narrower, BOLD impulse responses than the canonical HRF. This may reflect differences in vasculature (or even face-related neural differences across regions).

On the other hand, note that most voxels in
the above F-contrasts also show a positive loading on the canonical HRF (ie the
previous [1 0 0] T-contrast), as can be revealed by Inclusive (or Exclusive)
masking of the relevant contrasts. This is because the loadings on the
derivatives reflect deviations ABOUT the canonical form (via a first-order

One can also confirm this by going to various voxels in the above F-contrasts, pressing "plot", "contrast estimates" and selecting the "Can+Tem+Dis" F-contrast. The three bars indicate the loadings (and 90% confidence intervals) on the three different basis functions. Note that a positive estimate for the temporal derivative corresponds to an earlier response than the canonical (and negative for later), while a positive estimate for the dispersion derivative corresponds to a narrower (less dispersed) response (and negative for wider).

The contrast images for the difference between famous and nonfamous faces happen to correspond to numbers 39-74 in M1i, ie:

con_0039.img (canonical HRF, subject 1)

con_0040.img (canonical HRF, subject 2)

...

con_0051.img (temporal derivative, subject 1)

con_0052.img (temporal derivative, subject 2)

...

con_0063.img (dispersion derivative, subject 1)

con_0064.img (dispersion derivative, subject 2)

...

The model ("one-way ANOVA") is the same as for the main effect analysis in Eg2.1 above (so repeat the same steps, except selecting con_0039.img to con_0074.img instead).

Incidentally, the reason for not including a constant term (see above) is also why one cannot really examine the main effect and differential effect together within a single SPM design matrix consisting of 6 groups (ie 2 conditions and 3 basis functions). One might imagine that, this way, one could test contrasts like [1 0 0 1 0 0] (main effect) and [1 0 0 -1 0 0] (differential effect) in the same design matrix, ie without need to estimate separate SPMs for each effect (as in Eg2.1 and Eg2.2 here). However, one would also want to model subject effects, to remove between-subject variance, as is conventional in repeated-measures ANOVAs (such subject effects are already subtracted out in the current analysis by virtue of the fact that the 1st-level contrasts are already differences between conditions for each subject). But including a further 12 columns for the 12 subject effects would again make the design matrix rank deficient, so requiring contrasts that sum to zero and preventing tests of the main effect.

In the contrast manager press
'Define new contrast' (select F).

Enter ["eye(3)"] (which in matlab
evaluates to [1 0 0; 0 1 0; 0 0 1]) as the contrast weights, and
"Can+Tem+Dis: Famous vs Nonfamous" as the name.

This contrast will capture any voxels that show some difference in the shape of the event-related responses to Famous and Nonfamous faces (at least shape differences that can be captured by the three basis functions).

At a FWE corrected p-value of 0.05, this contrast reveals mainly a medial parietal cluster (3 -48 12, Z=5.49).

It is also informative to evaluate
the T-contrast [1 0 0] (ie positive loadings on the difference in canonical
HRFs, ie "Famous > Nonfamous").

Like with the more general F-contrast above, at p<0.05 corrected, this contrast reveals mainly a medial parietal cluster (3 -48 12, Z=6.12), though the higher Z-value indicates that the difference in response to famous and nonfamous faces in this region reflects mainly a difference in magnitude (of the canonical HRF).

The SPM is quite different from the equivalent contrast ([1]) in the model that uses only the canonical HRF (in Eg1.2 above). The main reason for this difference again relates to the different degrees of freedom: in the "canonical HRF only" model in Eg1.2, the 11 df's provide poor estimators of the underlying error, producing "speckled" SPMs (see above).

Again further information can be obtained by going to various voxels in the above contrasts, pressing "plot", "contrast estimates" and selecting the "Can+Tem+Dis Famous vs Nonfamous" F-contrast. The three bars indicate how the responses to the two conditions differ (loosely speaking, in terms of "magnitude", "latency" and/or "dispersion" of a canonical HRF - in this case, mainly in magnitude).

For this example, 12 contrast images per subject are taken to the 2nd-level.

For the main effect versus baseline, these are the contrast images:

con_fir_bin01_sub01.img (FIR bin 1, subject 1)

con_fir_bin01_sub02.img (FIR bin 1, subject 2)

...

con_fir_bin02_sub01.img (FIR bin 2, subject 1)

...

(available here: cons_fir.zip)

These images comprise the data for M2f, which is a "one-way ANOVA" in SPM2.

Type SPM at the matlab prompt.

Now change to a new directory (this is easy
to forget !)

Press the 'Basic models' button.

Select design type ..... [One way ANOVA]

#group's ... enter [12]

group 1: change the filter to
"con_fir_bin01*.img" and select "all"

group 2: change the filter to
"con_fir_bin02*.img" and select "all"

etc

GMsca: grand mean scaling [No]

Threshold masking [none]

explicitly mask images [No]

Global calculation [omit]

non-sphericity correction [yes]

replications are over?... [repl(12)]

correlated repeated measures [yes]

The "repl(12)" answer to the "replications are over?" question tells SPM2 that the replications correspond to subjects (SPMs implicit factor "replications"), with the other factor "group" corresponding to the repeated measure, namely the FIR timebins.

As before you can examine SPM.xVi.Vi

After estimation, test the "effects of interest" F-contrast: (Again, note the effect of pre-whitening on the design matrix)

In the contrast manager press 'Define
new contrast' (select F).

Enter ["eye(12)"] (which in matlab evaluates to a 12x12 identity matrix) as the contrast weights, and "FIR" as the name.

This contrast will reveal voxels that show ANY form (shape) of event-related response (within the range 0-24s poststimulus, and with 2s "resolution") (e.g, 42 -57 -24, Z=Inf (beyond SPM's range)).

Selecting a voxel and plotting this contrast will reveal that most voxels have a fairly "canonical" shape over the 12 timebins.

One can also test for more constrained shapes of event-related responses within this model. For example, one can test for "canonical-shaped" responses by evaluating a contrast whose weights trace out SPM's canonical HRF (every 2s).

To do this, switch to the Matlab window for a moment and type:

xBF.dt = 1;

xBF.name = 'hrf (with time and dispersion
derivatives)';

xBF.length = 32;

xBF.order = 1;

xBF = spm_get_bf(xBF);

This returns the canonical and two derivatives in the matrix "xBF.bf" (type "help spm_get_bf" for more info). For convenience, define:

can = xBF.bf(2:2:24,1)';

tem = xBF.bf(2:2:24,2)';

dis = xBF.bf(2:2:24,3)';

all = [can; tem; dis];

If you type "corrcoef(all')", you will see that the basis functions are slightly correlated (in the off-diagonal terms), due to their undersampling every 2s.

In the contrast manager press 'Define
new contrast' (select T).

Enter ["can"] as the contrast
weights (defined in Matlab workspace

as above), and "Can-weighted FIR" as the name.

At a FWE correct p value of 0.05, note more voxels (and higher Z-values e.g, 42 -57 -27, Z=Inf) for this main effect vs baseline compared to the equivalent contrasts in examples Eg1.1 and Eg2.1 above. The main reason for this increased power is again the increase in the degrees of freedom, which entails better estimators of the underlying error (co)variance (though if the FIR parameters were estimated very inefficiently, the extra contrast images might add more noise, outweighing any advantage of higher degrees of freedom). Again, this increased power comes with a stronger assumption about the nonsphericity, namely that it has the same structure across (activated) voxels – the pooling device" (see Glaser et al. (2003) Ch9.pdf).

One can also test the variance captured by the temporal and dispersion derivatives by creating new contrasts (though perhaps better as an F rather than T contrast) and simply typing "tem" and "dis" respectively as the contrast weights.

More interesting is the ability to ask, within this model, how much event-related variance is NOT captured by the canonical HRF. To do this, first create the variable in Matlab:

nullcan = eye(12) - pinv(can)*can;

This creates a matrix for an F-contrast that spans the "null space" of the canonical HRF.

In the contrast manager press 'Define
new contrast' (select F).

Enter ["nullcan"] as the contrast
weights (defined in Matlab workspace

as above), and "Null space of
canonical HRF" as the name.

You will see that several regions express variability not captured by the canonical HRF. This is not surprising, because you will notice that many of these regions appeared in the individual F-tests on the temporal and dispersion derivatives above, suggesting that what is not captured by the canonical HRF is captured by its two derivatives.

So even more interesting is the ability to ask how much event-related variance is NOT captured by the canonical HRF or its two derivatives (ie not captured by SPM's "informed" basis set). To do this, first create the variable in Matlab:

nullall = eye(12) - pinv(all)*all;

This creates a matrix for an F-contrast that spans the "null space" of all three informed basis functions.

In the contrast manager press 'Define
new contrast' (select F).

Enter ["nullall"] as the contrast
weights (defined in Matlab workspace

as above), and "Null space of informed
basis set" as the name.

You will see that only 2 voxels (in one cluster with maximum -21 -18 27) express variability not captured by the informed basis set. This reinforces the point that, while there is certainly variability in the HRF across different brain regions, the canonical HRF and its two derivatives are sufficient to capture the majority of this regional variability (at least on average across the 12 subjects in this dataset). See Henson et al (2000) hbm-fir.pdf for further details.

For the differential effect between famous and nonfamous faces, the contrast images are:

con_fir_dif_bin01_sub0001.img (FIR bin 1, subject 1)

con_fir_dif_bin01_sub0002.img (FIR bin 1, subject 2)

...

con_fir_dif_bin02_sub0001.img (FIR bin 2, subject 1)

...

These can be entered into a "one-way ANOVA" with 12 groups (like Eg3.1 above).

Here, evaluating the (effects of interest) F-contrast:

In the contrast manager press 'Define
new contrast' (select F).

Enter ["eye(12)"] (which in
matlab evaluates to a 12x12 identity matrix) as
the contrast weights, and "FIR Famous vs Nonfamous" as the name.

produces similar results to the equivalent differential contrast using only the three informed basis functions (in Eg2.2): The only cluster is in medial parietal cortex (0 -48 15, Z=5.73), which is bigger than the equivalent cluster in Eg2.2, but there are no other voxels, unlike the few other regions in Eg2.2. This again suggests that, when responses to famous and nonfamous faces differ, they differ mainly in a space captured by the informed basis set (ie not in other strange ways).

One can also test for differences between famous and nonfamous faces that are constrained to the canonical HRF (ie magnitude differences, loosely speaking). Assuming that you still have the variable “can” in Matlab workspace that was defined in Eg3.1 above:

In the contrast manager press 'Define
new contrast' (select T).

Enter ["can"], and "FIR
Canonical-weighted Famous > Nonfamous" as name.

The medial parietal cluster is now bigger and has a greater maximal value (0 -48 15, Z=6.71), illustrating that the difference between famous and nonfamous faces in this region is mainly in the magnitude of a canonical-like response.