SensorSpm - Meg Wiki

Revision 35 as of 2009-06-26 13:07:27

Clear message
location: SensorSpm

3D Sensor SPMs

3D Topography x time sensor SPMs can help you identify the latency and location of effects in sensor space. For example, here's a word-pseudoword effect (N=19) on magnetometers and gradiometers (RMS), thresholded for height (p<.001 unc) and extent (p<.05 corr):

attachment:WdNw3DSensorSpm.jpg

Creating 3D Sensor SPMs

1. Select an averaged file (m*.mat), split it into separate files for magnetometers, gradiometers, and gradiometer-magnitude (m*-mags.mat, m*-grds.mat, and m*-grms.mat, respectively). (Note you can do a similar analysis within a single subject by writing epochs to images).

clear S
S.D='mae_mydata.mat';
spm_eeg_splitFIF_grms(S);

2. Write image volumes.

% Select options (see help for spm_eeg_convertmat2ana3D):
clear S
S.interpolate_bad = 0;
S.n = 32;
S.pixsize = 3;
% Select trial types:
S.trialtypes = [1 2];
% Mags:
S.Fname = 'mae_mydata-mags.mat';
spm_eeg_convertmat2ana3D(S);
% Grad magnitude:
S.Fname = 'mae_mydata-grms.mat';
spm_eeg_convertmat2ana3D(S);

3. Smooth image volumes (optional, but good idea for grms, since distribution of RMS of gradiometer values [over voxels] is not Gaussian, even if by Central Limit Theorem, the distribution of the mean over enough subjects will be).

% Select files (or pass from stage above)
P = spm_select(Inf,'image','Select image file(s)',[],pwd,'.*img');
% Smoothness in x (mm), y (mm) and z (ms)
SmoothKernel = [5 5 10];
for n=1:size(P,1);
           [pth,nam,ext] = fileparts(P(n,:));
           Pout{n}       = fullfile(pth,['s' nam ext]);
           spm_smooth(spm_vol(P(n,:)),Pout{n},SmoothKernel);
           Pin = strvcat(P(n,:),Pout{n});
           spm_imcalc_ui(Pin,Pout{n},'((i1+eps).*i2)./(i1+eps)',{[],[],'float32',0});
end
%The final imcalc step above is just to reinsert NaN's for voxels outside space-time volume into which data were smoothed

4. Then compute an ANOVA over subjects. Either use SPM's GUI, or a script like 'meg_batch_anova.m' which can be found in /imaging/local/meg_misc (or here http://imaging.mrc-cbu.cam.ac.uk/svn/meg_misc/devel/), eg:

% wd is your working directory (eg pwd)
outdir = fullfile(wd,'SensorTimeSPMs');
meg_batch_anova(Pout,outdir);

Viewing 3D Sensor SPMs

To view the results, navigate to the directory that contains the SPM.mat for your ANOVA. In an SPM window, click 'Results' (1), select SPM.mat (2), select contrast Unwhitened EOI (3), and select options: (4) don't mask, whatever title, choose your p-value and extent thresholds (I usually start with .001 and 0).

attachment:Demo3dSpm1Scaled.jpg

Click 'whole brain' (5) to view whole-volume results. The coordinates on the far right (6) show x (mm), y (mm), and time (ms). Ignore the glass brain images (7).

attachment:Demo3dSpm2Scaled.jpg

To view the cloud-in-a-box images, click 'overlays'->'sections', (8) and select mask.img (9). Lovely (10). If you want to write out the image volume, click 'save' and give it a name. You can then view it using 'Display' or 'Check Reg'.

attachment:Demo3dSpm3Scaled.jpg

If you want to correct for extent, rather than clicking 'Results', click 'Toolbox'->'ns'. As above, select SPM.mat, Unwhitened EOI, options: don't mask, whatever title, no p-value adjustment, threshold {p} .001 (or your choice), extent threshold 0. The leftmost column ('cluster-level p-corrected') shows p-values corrected for cluster extent (11). Note that this information was absent when you clicked 'Results'.

attachment:Demo3dSpm4Scaled.jpg

If you want to view or write out an image thresholded for cluster extent, determine what cluster size ('cluster-level k-voxel') needed to achieve your desired p-value ('cluster-level p-corrected'; in the example, 214 gets p=0.052 (12)). Then follow the above 'Toolbox'->'ns' directions, but this time enter the k-voxel extent threshold when prompted (eg., 215 (13)). Now the results are thresholded for both height and extent (14). View and save as above.

attachment:Demo3dSpm5Scaled.jpg

A priori timewindows/sites of interest

Finally, if you have an a priori timewindow and/or an a priori set of sensors where you expect an effect, you can use Small Volume Correction (the SVC button). With this, you can specify a mask image saved from an orthogonal contrast in the same SPM analysis, or from a previous SPM analysis of independent data. Or perhaps more commonly, you can select a 3D "box". For example, say you had an a priori hypothesis that a condition effect would occur from 100 to 200ms poststimulus, then type [0 0 150] into the coordinates on the Input window, press SVC, select "box" and enter [80 80 100] as its dimensions (corresponding to an +/-40mm range around x, +/-40mm range around y, and +/-50ms range around 150ms). The corrected (but not uncorrected) p-values in the Graphics window will now change (get smaller) owing to the reduced search volume. Hopefully some will be <.05 corrected for height or extent!

A note on Gradiometer RMS

Sensor-level SPMs apply the same logic ([http://imaging.mrc-cbu.cam.ac.uk/imaging/PrinciplesRandomFields?highlight=%28field%29%7C%28random%29%7C%28theory%29 random-field theory, RFT]) as for fMRI data to MEG data, except that the image "volumes" do not consist of "voxels" in 3D-space, but rather have a 2D spatial dimension (locations of sensors) plus time as a 3rd dimension. For magnetometers, where there is one measurement per location, RFT can be applied directly. For gradiometers, however, there are two values at each location, because the gradiometer data reflect in-plane gradients (x,y) in two orthogonal directions. It is not straightforward to apply RFT (and e.g. define criteria for "smoothness") to such "vector fields". One could perform SPMs on each gradient direction separately, but this direction varies across locations with respect to the brain, so is not particularly meaningful. Furthermore, if the signal were present in both x and y directions, separate tests may be less sensitive than a combined test (and we would have to deal with potentially contradicting results from 3 sensor types, i.e. magnetometers and the two gradiometer types). It is therefore desirable to find a way of combing the two gradiometers values at each location. A common way of doing this the "Root-Mean_Square" (RMS): rms(x) = sqrt(sum(x₁²)/N), or for each gradiometer simply rms(g) = sqrt(sum(g²+g²)/2). This "rectification" of the signal removes information about its sign, or the direction of the gradient at each location. This is illustrated in the following images, where the left one respresents the original measurements of a gradiometer array (each arrow representing the two values of two gradiometers at one location), while the right one shows the RMS distribution computed from this vector field. The RMS values only reflect the length of the arrows, but not their orientation.

attachment:VectorFields.jpg

However, RMS (or any other rectification procedures) skews the distribution of original values toward positive values. Therefore, noise will not cancel when averaging the RMS (roughly speaking, positive and negative values do not cancel each other, but actually add up). To see this, consider that each measurement, y, is a mixture of signal (s) and noise (n), i.e, y=s+n (under additive noise assumptions). Let's assume, as typical, that the noise is zero-mean and Gaussian. Therefore, when averaging across many trials or subjects (observations of y), the expectation value <n>=0, and <y> approaches the true signal, s (in fact, the Signal-to-Noise ratio, SNR, increases with the square root of the number of trials). However, when you compute the signal power (part of the RMS), then the noise no longer cancels across trials (more formally, y² = (s+n)² = s² + s×n + n², and though <n>=0 (and <s×n>=0), n² is not zero). This means that even if there were no signal, the distribution of y across trials or subjects would be non-central (i.e, mean greater than zero). Thus one cannot test the null hypothesis that there is no signal (e.g, tests of an evoked peak versus baseline, when original data have been baseline-corrected); likewise, the RMS during the baseline period will be always be greater than zero. One might suggest "re-baseline-correcting" after taking the RMS, but this still makes assumptions about stationary noise, as explained below.

Because the RMS of the noise does not tend to zero with averaging, when comparing the RMS of two observations, where y=s+n and y=s+n, one has to assume that the noises, n and n, are equivalent. When comparing the same peristimulus timepoints across two trial-types, and assuming those trial-types are not confounded by some other variable such as time during experiment (unlikely if randomly intermixed), then the assumption that n=n would seem reasonable. However, when comparing two different peristimulus timepoints within the same trial, e.g, a pre-stimulus (baseline) time (where it is assumed s=0 and y=n) vs a post-stimulus time (y=s+n), the assumption that n=n is arguably less tenable, because the noise is unlikley to be stationary, e.g, because of random drifts through the trial.

Even if you think the assumption that n=n is valid, it is also important to remember that the distribution of RMS across any observations will be positively-skewed (non-Gaussian). In theory, this would invalidate parametric statistical tests (though some are quite robust, and the skewness depends on the ratio of fixed signal to random noise). So you might want to apply a further transform, e.g, log. One situation when this becomes more complex is when your stats are over subjects, but each subject has multiple trials. Normally one would take the RMS of the average across trials (for each subject), since the averaging will reduce the amount of noise (i.e, increase the SNR, as explained above). Therefore the distribution across subjects of this RMS value will be non-Gaussian (skewed), and a transform may be worthwhile. However, one could also calculate the RMS of each trial, and then average those within each subject. In this case, the noise will not cancel (SNR will be worse). On the positive side however, the distribution of the *mean* of any distribution, with enough observations, will approach a Gaussian (by the Central Limit Theorem). Therefore if your statistical tests are over enough subjects, then parametric tests on the mean-across-trials of the RMS will be valid.

For further information about the gradiometer RMS issue, see SensorStats.