Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revision Previous revision
Next revision
Previous revision
example:regressing_out_headposition_confounds_in_a_ctf275_dataset [2012/05/29 17:13]
131.174.45.254 [Description]
example:regressing_out_headposition_confounds_in_a_ctf275_dataset [2017/08/17 11:21] (current)
Line 1: Line 1:
-{{tag>​example artifact meg regression confound}} +Check here for new version ​of this topic [[example:how_to_incorporate_head_movements_in_meg_analysis|How to incorporate ​head movements ​in MEG analysis]]
-====== How to incorporate head movements in MEG analysis ====== +
- +
-===== Description ===== +
- +
-Changes in head position during MEG session may cause significant error in the source localization. Besides, the mixture ​of different head positions over time adds variance to the data that is not accounted for by the experimental manipulation. Thus head movements may deteriorate statistical sensitivity when analyzing MEG on both sensor and source levels. It is recommended to incorporate head movements in MEG analysis. +
- +
-Continuous head localization information is stored in HLC channels (Head Localization Channels) in CTF MEG system. An example script here shows how to read these channels in FieldTrip and estimate the amount of movement offline. The information from these channels can also be used to track the head position in real time [[faq:​how_can_i_monitor_a_subject_s_head_position_during_a_meg_session|(see here).]] +
-The data used in this example ​script ​ [[ftp://​ftp.fcdonders.nl/​pub/​fieldtrip/​tutorial/​TactileStimulusDipolefit.zip|can be obtained here]] +
- +
- +
- +
-===== Reading-in HLC channels ===== +
- +
-Prepare configuration ​to define trials:  +
- +
-<code matlab>​ +
- +
-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);​ +
- +
-</​code>​ +
- +
-Read the data with the following HLC channels:  +
- +
-  * HLC00n1 X coordinate relative to the dewar (in meters) of the nth head localization coil +
-  * HLC00n2 Y coordinate relative to the dewar (in meters) of the nth head localization coil +
-  * HLC00n3 Z coordinate relative to the dewar (in meters) of the nth head localization coil +
- +
-<code matlab>​ +
-cfg.channel ​                = {'​HLC0011','​HLC0012','​HLC0013',​ ... +
-                              '​HLC0021','​HLC0022','​HLC0023',​ ... +
-                              '​HLC0031','​HLC0032','​HLC0033'​};​ +
- +
-headpos = ft_preprocessing(cfg);​ +
-</​code> ​  +
- +
-Determine the mean (per trial) circumcenter (the center of the circumscribed circle) of the three headcoils and its orientation (see subfunction) +
- +
-<code matlab>​ +
-ntrials = length(headpos.sampleinfo);​ +
-for t = 1:ntrials +
-  % x,y,z of coil 1, 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 +
-</​code>​ +
- +
-Now you can plot the head position relative to the first value, and compute the maximal position change. +
- +
-<code matlab>​ +
- +
-cc_rel = [cc - repmat(cc(:,​1),​1,​size(cc,​2))]';​ +
- +
-% plot translations +
-figure();  +
-plot(cc_rel(:,​1:​3)*1000) % in mm +
- +
-% plot rotations +
-figure();  +
-plot(cc_rel(:,​4:​6)) +
- +
-maxposchange = max(abs(cc_rel(:,​1:​3)*1000)) % in mm +
- +
-</​code>​ +
- +
-You may decide to exclude a subject from the subsequent analysis if the head movement exceeds a certain threshold. You may also incorporate the movement estimate as a regressor in your analysis. +
- +
- +
-===== Regressing out headposition confounds ===== +
- +
-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.  +
- +
-<code matlab>​ +
-% 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); +
-</​code>​ +
- +
-Use the circumcenter,​ as computed in the paragraph above. +
- +
-<code matlab>​ +
-% 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);​ +
-</​code>​ +
-===== Figure ===== +
-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). +
- +
-{{:​example:​regr_ftwiki.png?&​300|}} +
- +
-===== Subfunction ===== +
-<code matlab>​ +
-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 +
-</​code>​+