Magnetoencephalography (MEG) measurements are strongly deteriorated by head movements, leading to reduced statistical sensitivity. Also, the mixture of different head positions over time adds variance to the data that is not accounted for by the experimental manipulation, thus potentially deteriorating statistical sensitivity. Offline incorporation of the head position time-series into the general linear model results in improvements of group-level statistical sensitivity.

Trial-by-trial estimates of the position of the circumcenter, i.e. the center of the circle that passes through all the coil positions, can be used as a regressor. These are then demeaned and z-transformed to obtain the normalized deviants, i.e. translations, from the average head position. The estimated contributions to the (source reconstructed) signal amplitude and frequency power are then removed from the single-trial data at sensor- and source level.

The Matlab script is given for an example timelock analysis pipeline. Note that ft_regressconfound can be applied to both sensor- and sourcelevel timelock/freq data. Use it as the last step prior to ft_XXXstatistics.

The MEG dataset is available from ftp://ftp.fcdonders.nl/pub/fieldtrip/example/regressconfound and is named TacStimRegressConfound.ds

% trial definition
cfg                         = [];
cfg.dataset                 = 'TacStimRegressConfound.ds';
cfg.trialdef.eventtype      = 'UPPT001';
cfg.trialdef.eventvalue     = 4;
cfg.trialdef.prestim        = 0.2;
cfg.trialdef.poststim       = 0.3;
cfg.continuous              = 'yes';
cfg = ft_definetrial(cfg);
 
% preprocess the headposition
cfg.channel                 = {'HLC0011','HLC0012','HLC0013', ...
                              'HLC0021','HLC0022','HLC0023', ...
                              'HLC0031','HLC0032','HLC0033'};
headpos = ft_preprocessing(cfg);
 
% preprocess the data
cfg.channel                 = {'MEG'};
cfg.demean                  = 'yes'; 
cfg.baselinewindow          = [-0.2 0];
cfg.dftfilter               = 'yes'; % notch filter to filter out 50Hz
data = ft_preprocessing(cfg);
 
% timelock analysis
cfg                         = [];
cfg.keeptrials              = 'yes';
timelock = ft_timelockanalysis(cfg, data);
 
% determine the mean (per trial) circumcenter of the three headcoils and its orientation
% (see subfunction)
ntrials = length(headpos.sampleinfo);
for t = 1:ntrials
  % x,y,z of coil1, 2, and 3
  cc(:,t) = circumcenter( ...
  [mean(headpos.trial{1,t}(1,:)); mean(headpos.trial{1,t}(2,:)); mean(headpos.trial{1,t}(3,:))], ...
  [mean(headpos.trial{1,t}(4,:)); mean(headpos.trial{1,t}(5,:)); mean(headpos.trial{1,t}(6,:))], ...
  [mean(headpos.trial{1,t}(7,:)); mean(headpos.trial{1,t}(8,:)); mean(headpos.trial{1,t}(9,:))]);
end
 
% demean to obtain translations and rotations from the average head
% positions and orientation
cc_dem = [cc - repmat(mean(cc,2),1,size(cc,2))]';
 
% add these, their squares, cubes and all their derivatives to
% the regressorlist. also add the constant (at the end; column 37)
confound = [cc_dem cc_dem.^2 cc_dem.^3 ...
gradient(cc_dem')' gradient(cc_dem.^2')' gradient(cc_dem.^3')' ...
ones(size(cc_dem,1),1)];
 
% regress out headposition confounds
cfg = [];
cfg.confound = confound;
cfg.reject = [1:36]; % keeping the constant (nr 37)
regr = ft_regressconfound(cfg, timelock);

Example statistical results in a single-subject (baseline vs. task activity contrasts). With ft_regressconfound, sensor-level statistical sensitivity was increased after tactile stimulation (40-50 ms; note the more extreme t-scores in the upper panel). In a similar vein, source-level statistical sensitivity was increased after visual stimulation (0-500 ms; 65Hz; lower panel).

function [cc] = circumcenter(coil1,coil2,coil3)
 
% CIRCUMCENTER determines the position and orientation of the circumcenter
% of the three fiducial markers (MEG headposition coils). 
%
% Input: X,y,z-coordinates of the 3 coils [3 X N],[3 X N],[3 X N] where N
% is timesamples/trials.
%
% Output: X,y,z-coordinates of the circumcenter [1-3 X N], and the 
% orientations to the x,y,z-axes [4-6 X N].
%
% A. Stolk, 2012
 
% number of timesamples/trials
N = size(coil1,2);
 
%% x-, y-, and z-coordinates of the circumcenter
% use coordinates relative to point `a' of the triangle
xba = coil2(1,:) - coil1(1,:);
yba = coil2(2,:) - coil1(2,:);
zba = coil2(3,:) - coil1(3,:);
xca = coil3(1,:) - coil1(1,:);
yca = coil3(2,:) - coil1(2,:);
zca = coil3(3,:) - coil1(3,:);
 
% squares of lengths of the edges incident to `a'
balength = xba .* xba + yba .* yba + zba .* zba;
calength = xca .* xca + yca .* yca + zca .* zca;
 
% cross product of these edges
xcrossbc = yba .* zca - yca .* zba;
ycrossbc = zba .* xca - zca .* xba;
zcrossbc = xba .* yca - xca .* yba;
 
% calculate the denominator of the formulae
denominator = 0.5 ./ (xcrossbc .* xcrossbc + ycrossbc .* ycrossbc + zcrossbc .* zcrossbc);
 
% calculate offset (from `a') of circumcenter
xcirca = ((balength .* yca - calength .* yba) .* zcrossbc - (balength .* zca - calength .* zba) .* ycrossbc) .* denominator;
ycirca = ((balength .* zca - calength .* zba) .* xcrossbc - (balength .* xca - calength .* xba) .* zcrossbc) .* denominator;
zcirca = ((balength .* xca - calength .* xba) .* ycrossbc - (balength .* yca - calength .* yba) .* xcrossbc) .* denominator;
 
cc(1,:) = xcirca + coil1(1,:);
cc(2,:) = ycirca + coil1(2,:);
cc(3,:) = zcirca + coil1(3,:);
 
%% orientation of the circumcenter with respect to the x-, y-, and z-axis
% coordinates
v = [cc(1,:)', cc(2,:)', cc(3,:)'];
vx = [zeros(1,N)', cc(2,:)', cc(3,:)']; % on the x-axis
vy = [cc(1,:)', zeros(1,N)', cc(3,:)']; % on the y-axis
vz = [cc(1,:)', cc(2,:)', zeros(1,N)']; % on the z-axis
 
for j = 1:N
  % find the angles of two vectors opposing the axes
  thetax(j) = acos(dot(v(j,:),vx(j,:))/(norm(v(j,:))*norm(vx(j,:))));
  thetay(j) = acos(dot(v(j,:),vy(j,:))/(norm(v(j,:))*norm(vy(j,:))));
  thetaz(j) = acos(dot(v(j,:),vz(j,:))/(norm(v(j,:))*norm(vz(j,:))));
 
  % convert to degrees
  cc(4,j) = (thetax(j) * (180/pi));
  cc(5,j) = (thetay(j) * (180/pi));
  cc(6,j) = (thetaz(j) * (180/pi));
end