Statistics in sensor space
The following scripts are based on SensorSpm, which offers a standardised way of running sensor-level statistics on your MEG data in SPM5. The following scripts may help you
- convert fiff-files to SPM images (if necessary)
- combine gradiometers such that condition can be tested against zero
- run group statistics on MEG data in signal space in SPM5
1. Converting Fiff-files to SPM images (for between-condition contrasts)
This script is based http://imaging.mrc-cbu.cam.ac.uk/meg/SensorSPM. It converts magnetometer and gradiometer data to SPM image files, and computes RMS values for gradiometers (2 at each location) if requested. This is appropriate for statistical contrasts between conditions, but not for comparisons of one condition against zero. If that's what you're after, look under point 2 below.
% Modified from http://imaging.mrc-cbu.cam.ac.uk/meg/SensorSPM
% takes a list of Fiff-files (cell array "fifflist{}") and converts MEG data to SPM image files,
% which can then be subjected to 2nd-level statistics, for example
% "rootdir" can contain the directory path common to all files
% the script will produce subdirectories with img-files in the directories of the input fiff-files
% so far only tested on one data set
% OH, March 2009
% you can specify a list of fiff-files here, to be converted into SPM format
fifflist = {'/fullpath/file4subj1.fif', ...
'/fullpath/file4subj2.fif', ...
'/fullpath/file4subj3.fif'};
nr_fiffs = length(fifflist);
% root directory common to all input fiff-files (not required if full pathnames specified in fifflist{})
root_dir = '';
% GET GOING...
clear S;
for ff = 1:nr_fiffs, % for all Fiff-files...
% CONVERT FROM FIFF TO SPM FORMAT (for between-condition contrasts)...
S.Fchannels = fullfile(spm('dir'),'EEGtemplates','FIF306_setup.mat');
S.Fchannels_eeg = fullfile(spm('dir'),'EEGtemplates','cbu_meg_70eeg_ecg_montage.mat');
S.veogchan = [62]; % Channel 62 is (bipolar) VEOG, MAY NEED EDITING
S.heogchan = [61]; % Channel 61 is (bipolar) HEOG, MAY NEED EDITING
S.ecgchan = [63]; % MAY NEED EDITING (e.g. you may not use ECG)
S.veog_unipolar = 0; % MAY NEED EDITING
S.heog_unipolar = 0; % MAY NEED EDITING
S.ecg_unipolar = 0; % MAY NEED EDITING (e.g. you may not use ECG)
S.dig_method = 1; % new digitising conventions
S.conds = 1; % conditions in fiff-file, MAY NEED EDITING
S.grms = 1; % so that the splitting outputs mags, grads and grad RMS (grms)
[fiffpath,fiffname,fiffext,fiffversn] = fileparts(fifflist{ff});
S.Fdata = fullfile(root_dir, fifflist{ff}); % Input Fiff-file before splitting
fileSPMout = fullfile(root_dir, fiffpath, [fiffname '_SPM.mat']); % Output SPM file (if not specified, will be 'myexp.mat')
S.Pout = fileSPMout; % output for conversion from Fiff to SPM format (before splitting)
D0 = spm_eeg_rdata_FIF(S); % convert fiff- to SPM-format
% SPLIT DATA INTO MAGS/GRADS/EEG...
S.D = fileSPMout; % structure for splitting
D1 = spm_eeg_splitFIF(S); % split SPM files
% WRITE TO 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];
img_outnames = '';
% Mags:
S.Fname = fullfile(root_dir, fiffpath, [fiffname '_SPM-mags.mat']);
tmpname = spm_eeg_convertmat2ana3D(S);
img_outnames(1,:) = fullfile(root_dir, fiffpath, [fiffname '_SPM-mags'], 'trialtype1', 'average.img');
% GradRMS magnitude:
S.Fname = fullfile(root_dir, fiffpath, [fiffname '_SPM-grms.mat']);
tmpname = spm_eeg_convertmat2ana3D(S);
img_outnames(3,:) = fullfile(root_dir, fiffpath, [fiffname '_SPM-grms'], 'trialtype1', 'average.img');
% EEG magnitude:
S.Fname = fullfile(root_dir, fiffpath, [fiffname '_SPM-eeg.mat']);
tmpname = spm_eeg_convertmat2ana3D(S);
img_outnames{3} = fullfile(root_dir, fiffpath, [fiffname '_SPM-eeg'], 'trialtype1', 'average.img');
% SMOOTH IMAGES
P = img_outnames;
SmoothKernel = [5 5 10]; % % Smoothness in x (mm), y (mm) and z (ms), MAY NEED EDITING
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
fnames_for_stats{ff} = Pout;
end; %..ff
2. Converting Fiff-files to SPM images (for single-condition contrasts against zero)
If you are testing a single condition against zero, computing the RMS for the two gradiometers at each sensor location increases the likelihood of false positives. Even the baseline noise will have RMS values above zero. One way of minimising this problem is computing the signal power for each gradiometer pair, and subtracting the mean baseline value afterwards. Unfortunately, this requires a different script at the moment:
% takes a list of Fiff-files (cell array "fifffiles{}") and converts MEG data to SPM image files
% different output options for gradiometers can be specified
% OH, June 2009
addpath /imaging/olaf/MEG/SPM_tools; % for file meg_grds2grms_pow.m, which computes power and baseline correction for grads
% you can specify a list of fiff-files here, to be converted into SPM format
fifffiles = {'/fullpath/file4subj1.fif', ...
'/fullpath/file4subj2.fif', ...
'/fullpath/file4subj3.fif'};
root_dir = ''; % path relative to which filenames are specified
%% parameters for 1) CONVERT FROM FIFF TO SPM FORMAT... (spm_eeg_splitFIF)
clear S;
S.Fchannels = fullfile(spm('dir'),'EEGtemplates','FIF306_setup.mat');
S.Fchannels_eeg = fullfile(spm('dir'),'EEGtemplates','cbu_meg_70eeg_ecg_montage.mat');
S.veogchan = [62]; % Channel 62 is (bipolar) VEOG, MAY NEED EDITING
S.heogchan = [61]; % Channel 61 is (bipolar) HEOG, MAY NEED EDITING
S.ecgchan = [63]; % MAY NEED EDITING (e.g. you may not use ECG)
S.veog_unipolar = 0; % MAY NEED EDITING
S.heog_unipolar = 0; % MAY NEED EDITING
S.ecg_unipolar = 0; % MAY NEED EDITING (e.g. you may not use ECG)
S.dig_method = 1; % new digitising conventions
S.conds = 1; % conditions in fiff-file, MAY NEED EDITING
S.grms = 1; % so that the splitting outputs mags, grads and grad RMS (grms)
%% parameters for 2) Compute RMS for grads, with/without log (meg_grds2grms_pow)
clear SR;
SR.do_write = 1; % write result files
SR.pow = 1; % compute power for grads
SR.log = 0; % take logarithm of grads
SR.bc = 1; % correct baseline after RMS/log
%% parameters for 3) convert fiff-files to SPM images (incl. smoothing) (spm_eeg_convertmat2ana3D, spm_smooth)
clear SC;
SC.interpolate_bad = 0;
SC.n = 32;
SC.pixsize = 3;
SC.trialtypes = [1];
%% MAIN PROCESSING...
fprintf(1, 'We are starting...\n\n');
ini_dir = pwd; % will go back to current working directory after processing
nr_fiffs = length(fifffiles);
%% next section processes both magnetometers and gradiometers
for ff = 1:nr_fiffs, % for all Fiff-files...
%% 1) CONVERT FROM FIFF TO SPM FORMAT...
[fiffpath,fiffname,fiffext,fiffversn] = fileparts(fifffiles{ff});
S.Fdata = fullfile(root_dir, fifffiles{ff}); % Input Fiff-file before splitting
S.Pout = fullfile(root_dir, fiffpath, [fiffname '_SPM.mat']); % Output SPM file (if not specified, will be 'myexp.mat')
fprintf(1, 'spm_eeg_rdata_FIF\n');
spm_eeg_rdata_FIF(S); % convert fiff- to SPM-format
%% SPLIT DATA INTO MAGS/GRADS/EEG...
S.D = S.Pout; % structure for splitting
fprintf(1, 'spm_eeg_splitFIF\n');
spm_eeg_splitFIF(S); % split SPM files
%% 2) Compute RMS for grads, with/without log
grad_file = fullfile( root_dir, fiffpath, [fiffname '_SPM-grds.mat'] ); % gradiometer files only, for later
SR.D = grad_file;
fprintf(1, '\nmeg_grds2grms_pow\n');
meg_grds2grms_pow(SR); % SPM script that splits data into mags and grads, computes power and baseline correction
%% 3) convert fiff-files to SPM images
fprintf(1, '\nWrite to image volumes\n');
clear P;
% Select options (see help for spm_eeg_convertmat2ana3D):
img_outnames = '';
% figuring out filenames for gradiomter data after conversion, depending on parameters
try dopow=SR.pow; catch dopow=0; end;
try dolog=SR.log; catch dolog=0; end;
try bc=SR.bc; catch bc=0; end;
namepart = '';
if dopow,
namepart = [namepart 'p'];
else,
namepart = [namepart 'r'];
end;
if dolog,
namepart = [namepart 'l'];
end;
if bc,
namepart = [namepart 'b'];
end;
switch length(namepart)
case 1,
namepart = ['-grd' namepart];
case 2,
namepart = ['-gr' namepart];
case 3,
namepart = ['-g' namepart];
end;
% Grad power magnitude:
SC.Fname = fullfile( root_dir, fiffpath, [fiffname '_SPM' namepart '.mat'] );
fprintf(1, 'spm_eeg_convertmat2ana3D\n');
tmpname = spm_eeg_convertmat2ana3D(SC);
P(1,:) = fullfile(root_dir, fiffpath, [fiffname '_SPM' namepart], 'trialtype1', 'average.img');
% Magnetometer magnitude
SC.Fname = fullfile( root_dir, fiffpath, [fiffname '_SPM-mags.mat'] );
fprintf(1, 'spm_eeg_convertmat2ana3D\n\n');
tmpname = spm_eeg_convertmat2ana3D(SC);
P(2,:) = fullfile(root_dir, fiffpath, [fiffname '_SPM-mags'], 'trialtype1', 'average.img');
% EEG magnitude:
SC.Fname = fullfile(root_dir, fiffpath, [fiffname '_SPM-eeg.mat']);
tmpname = spm_eeg_convertmat2ana3D(SC);
P(3,:) = fullfile(root_dir, fiffpath, [fiffname '_SPM-eeg'], 'trialtype1', 'average.img');
%% SMOOTH IMAGES (necessary for random-field theory to apply - still a bit of a mystery)
SmoothKernel = [5 5 10]; % % Smoothness in x (mm), y (mm) and z (ms), MAY NEED EDITING
for n=1:size(P,1);
[pth,nam,ext] = fileparts(P(n,:));
Pout{n} = fullfile(pth,['s' nam ext]);
disp('spm_smooth');
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});
disp(Pout{n});
end
fnames_for_stats{ff} = Pout;
%The final imcalc step above is just to reinsert NaN's for voxels outside space-time volume into which data were smoothed
end; %..ff
cd( ini_dir );
fprintf(1, '\nDone enough.\n');
3. Sensor-level group analysis in SPM 5
% based on http://imaging.mrc-cbu.cam.ac.uk/meg/SensorSpm
% to be run in SPM 5 EEG
% runs group statistics on files in "imgfiles"
% output written into directory "outdir"
% OH, March 2009
addpath /imaging/local/meg_misc; % for meg_batch_anova
% output directory
outdir = '/youroutputwillgohere/;
% Here: One group of subjects, 2 conditions
% You can also specify more groups of subjects, or test one condition against zero:
% More groups of subjects: Vary first cell index of imgfiles (imgfiles{1}{cc}...imgfiles{2}{cc}...)
% One-sample t-test against zero mean: just specify one file per line,
% e.g. imgfiles{1}{cc} = ['/thispath/subj1_con1.img']; cc=cc+1; etc.
if exist('imgfiles')~=1,
cc = 1;
imgfiles{1}{cc} = ['/thispath/subj1_con1.img'; 'thispath/subj1_con2.img']; cc=cc+1;
imgfiles{1}{cc} = ['/thispath/subj2_con1.img'; 'thispath/subj2_con2.img']; cc=cc+1;
imgfiles{1}{cc} = ['/thispath/subj3_con1.img'; 'thispath/subj3_con2.img'];
end;
nr_subjects = cc;
% root directory common to all input fiff-files (not required if full pathnames specified in imgfiles{})
root_dir = '';
% Attach root directory
clear Panova;
for i=1:nr_subjects, % create full image file names, including root directory etc.
[m,n] = size(imgfiles{1}{i});
for j=1:m,
Panova{1}{i}(j,:) = fullfile( root_dir, imgfiles{1}{i}(j,:) );
end; %..j
end; %..i
meg_batch_anova(Panova,imgset(ss).outdir); % Run SPM stats