[FieldTrip] [Eeglablist] nested hypothesis testing to decide whether to use one or two dipoles to fit a component
Maximilien Chaumon
maximilien.chaumon at gmail.com
Fri Aug 24 17:13:22 CEST 2012
Hi everyone,
Just to follow-up on this if anyone wanted to take that approach.
I did use the following procedure to decide whether to fit one or two
dipoles to a given component.
1) fit one dipole, compute SSE for each component map as the RV * sum of
squared inverse weights (sse(1) = EEG.dipfit.model(idip).rv .*
sum(EEG.icawinv(:,idip).^2))
2) fit two dipoles, compute SSE for each component map: sse(2)
3) assume p1 = 6 parameters estimated for 1 dipole, p2 = 9 parameters
estimated for 2 dipoles (only three momentum values more, since symmetry
imposes position for the two dipoles), n = EEG.nbchan
4) F statistic : F = (sse(1)-sse(2))/(p2-p1)/(sse(2)/(n-p2))
5) if F > finv(.95,p2-p1,n-p2) decide to go for 2 dipoles.
6) that yielded some solutions where 2 dipoles were selected but were in
fact at the same position (x == 0) with opposite momenta. Some others had
really just one dipole with large amplitude and the other much smaller. So
I discarded solutions where distance between the two dipoles was less than
5mm, or where the momentum angles were opposed, or when the amplitude ratio
between the two dipoles was bigger than 5.
The final solution is only moderately satisfying. Looking at the maps, the
2 dipole solution is still too often favored over the one dipole solution.
As Scott previously mentioned, it could be that the second dipole accounts
for error in the head model, and that's all. The F difference in SSE is
probably not a good measure of the two-dipolarity of a component.
One would probably better just screen all components by eye and spot those
that look dipolar using thumbometric visual assessment.
Best,
Max
PS: the script I used:
dofit = 1;
dorework = 0;
doplot = 1;
rootdir = '/DATA/';
cd(rootdir)
studyfile = fullfile(rootdir,'AP2.study');
if dofit || dorework
[STUDY ALLEEG] = pop_loadstudy(studyfile);
for i_set = 1:numel(STUDY.datasetinfo)
EEG =
pop_loadset('filename',ALLEEG(i_set).filename,'filepath',ALLEEG(i_set).filepath);%pop_loadset(GRAB(i_set).name);
if isfield(EEG,'peakalphafrequency')
EEG.peakalphafreq = EEG.peakalphafrequency;
EEG = rmfield(EEG,'peakalphafrequency');
end
if ~dorework
EEG = pop_dipfit_settings( EEG, 'hdmfile',[
fileparts(which('eeglab'))
'/plugins/dipfit2.2/standard_BEM/standard_vol.mat'],'coordformat','MNI','mrifile',[
fileparts(which('eeglab'))
'/plugins/dipfit2.2/standard_BEM/standard_mri.mat'],'chanfile',[
fileparts(which('eeglab'))
'/plugins/dipfit2.2/standard_BEM/elec/standard_1005.elc'],'coord_transform',[0.67468
-15.6307 2.426 0.081417 0.0024349 -1.5728 1.1744 1.0622 1.1487]
,'chansel',[1:EEG.nbchan] );
EEGonedipole = pop_multifit(EEG, [1:size(EEG.icawinv,2)]
,'threshold',100,'rmout','on','dipoles',1);
EEGtwodipoles = pop_multifit(EEG, [1:size(EEG.icawinv,2)]
,'threshold',100,'rmout','on','dipoles',2);
else
EEGonedipole = EEG;
EEGtwodipoles = EEG;
EEGonedipole.dipfit.model = EEG.dipfit.onedipole;
EEGtwodipoles.dipfit.model = EEG.dipfit.twodipoles;
end
% H0: one dipole.
EEG.dipfit = EEGonedipole.dipfit;
% still we store them both. doesn't hurt.
EEG.dipfit.onedipole = EEGonedipole.dipfit.model;
EEG.dipfit.twodipoles = EEGtwodipoles.dipfit.model;
% now for every component, if the rv is much better for the
% 2dipoles than for the 1dipole, we'll use the 2dipoles.
for i_dip = 1:numel(EEGonedipole.dipfit.model)
sse(1) = EEGonedipole.dipfit.model(i_dip).rv .*
sum(EEGonedipole.icawinv(:,i_dip).^2);
sse(2) = EEGtwodipoles.dipfit.model(i_dip).rv .*
sum(EEGtwodipoles.icawinv(:,i_dip).^2);
p1 = 6;% six params for 1dipole model
p2 = 9;% 9 for 2 dipoles model (assume only symmetry is fixed).
n = EEG.nbchan; % N is number of electrodes
% F statistic
F = (sse(1)-sse(2))/(p2-p1)/(sse(2)/(n-p2));
EEG.dipfit.Fone2two(i_dip) = F;
% F decision
Fthresh = finv(.999,p2-p1,n-p2);
if F > Fthresh % reject H0, take 2 dipoles model
% only if x ~= 0 and the two momenta are not just opposite.
mom = EEGtwodipoles.dipfit.model(i_dip).momxyz;
pos = EEGtwodipoles.dipfit.model(i_dip).posxyz;
[t p r] = cart2sph(mom(:,1),mom(:,2),mom(:,3));
t(2) = t(2)+pi;%
d = sqrt(sum(diff(pos).^2));
% if the two dipoles are too close and have exactly
% opposite angles
cond = d < 5 && (diff(abs(t)) < 0.01 || diff(abs(p)) <
0.01);
% or if they are exactly on the sagittal plane
cond = cond |
abs(EEGtwodipoles.dipfit.model(i_dip).posxyz(1)) < 5;
% or if the norm of one is much bigger than the other
ratio = norm(mom(1,:)) ./ norm(mom(2,:));
ratiothresh = 5;
cond = cond | (ratio > ratiothresh || ratio <
1/ratiothresh);
% we will still keep the one dipole fit
if ~cond% so if not, we take the 2 dipole fit.
EEG.dipfit.model(i_dip) =
EEGtwodipoles.dipfit.model(i_dip);
end
end
end
% reject poor overall fits.
EEG.dipfit.model = dipfit_reject(EEG.dipfit.model, 20/100);
% reject NaN rv
for i_m = 1:numel(EEG.dipfit.model)
if isnan(EEG.dipfit.model(i_m).rv)
EEG.dipfit.model(i_m).posxyz = [];
EEG.dipfit.model(i_m).momxyz = [];
EEG.dipfit.model(i_m).rv = 1;
end
end
% removing dipoles outside the head
% ---------------------------------
rmdip = [];
for index = 1:numel(EEGonedipole.dipfit.model)
if ~isempty(EEG.dipfit.model(index).posxyz)
if any(sqrt(sum(EEG.dipfit.model(index).posxyz.^2,2)) > 85)
rmdip = [ rmdip index];
EEG.dipfit.model(index).posxyz = [];
EEG.dipfit.model(index).momxyz = [];
EEG.dipfit.model(index).rv = 1;
end;
end;
end;
if ~isempty(rmdip)
fprintf('%d out of cortex dipoles removed (usually
artifacts)\n', length(rmdip));
end;
pop_saveset(EEG,'savemode','resave');
end
end
%%
if doplot
% [STUDY ALLEEG] = pop_loadstudy('AlphaPsyc2.study');
num = [];% number of dipoles that are fit for each component
idx = [];
topos = {[] [] []};
Fdiff = {[] [] []};
idx = {[] [] []};
dips = {};
for i_e = 1:numel(ALLEEG)
for i_comp = 1:numel(ALLEEG(i_e).dipfit.model)
num(end+1) = size(ALLEEG(i_e).dipfit.model(i_comp).posxyz,1);
topos{num(end)+1}(end+1,:) = ALLEEG(i_e).icawinv(:,i_comp)';
dips{num(end)+1}(size(topos{num(end)+1},1)) =
ALLEEG(i_e).dipfit.model(i_comp);
Fdiff{num(end)+1}(end+1) = ALLEEG(i_e).dipfit.Fone2two(i_comp);
idx{num(end)+1}(end+1,:) = [i_e i_comp];
end
end
for i_plot = 2:3
figure(6548+i_plot);paste_figpos manytopos
ntopos = size(topos{i_plot},1);
n = ceil(sqrt(ntopos));
m = n;
for i_topo = 1:ntopos
subplot(n,m,i_topo)
cla
topoplot(topos{i_plot}(i_topo,:),ALLEEG(1).chanlocs,'conv','on','electrodes','off');
title({num2str(idx{i_plot}(i_topo,:))
num2str(Fdiff{i_plot}(i_topo))});
drawnow
end
end
end
2012/8/13 Maximilien Chaumon <maximilien.chaumon at gmail.com>
> Hello eeglab & Fieldtrip,
>
> I'm trying to find out if it would be possible to use a nested hypothesis
> testing approach to decide whether to use a one or two dipole model while
> estimating components' dipole locations.
> The rationale I would like to follow is this: with two dipoles, we will
> always obtain a better fit than with one dipole, but the decrease in sum of
> squared errors (SSE) should follow a F distribution with k (=
> Nparameters_2dipoles - Nparameters_1dipole) degrees of freedom. If the
> decrease in SSE is greater than what would be expected under this F
> distribution, then we decide that 2 dipoles provide a sufficiently better
> fit and decide using them.
>
> I asked this question to eeglablist and Scott pointed out that it is
> difficult/impossible(?) to determine if the second dipole fits actual
> interesting data or just noise introduced by the imperfect head model.
> Christian then said it'd be worth a shot, and I agree, so here I am again
> with two questions, or two confirmations:
>
> 1) *How many parameters are estimated in ft_dipolefitting.m ?* specially
> in the case of 2 dipoles. If I count correctly, we estimate 6 parameters
> for one dipole, and, depending on whether the orientation has to be the
> same in the 2 dipoles, one (amplitude) or three (amplitude and orientation)
> more.
> 2) *Can I assimilate the relative residual variance to a SSE?* the
> function rv.m does this: rv = sum((d1-d2).^2) ./ sum(d1.^2);
> So that seems to be a sum of squared errors divided by the variance of the
> original data. So if I multiply the rv by the sum squared component map, I
> should get it, right?
>
> Thanks a lot!
> Max
>
>
> 2012/8/11 Christian Kothe <christiankothe at gmail.com>
>
>> I can only speak from my armchair here, but it sounds like it should be
>> worth a try - even if you don't get the # of parameters exactly right it
>> will probably give you at least some level of complexity control in
>> whatever the range of validity is. If it works, it may inspire follow-up
>> work (e.g., Bayesian model selection or likelihood ratio tests).
>>
>> The number of parameters for a 2-dipole model seems to be 3 (xyz) + 4
>> (2x the orientation parameters). Not sure about the momentum, though - you
>> might look up the place where the actual function minimization is being
>> performed in dipfit (fminunc call?) and see whether these are being
>> optimized together with the others.
>>
>> Christian
>>
>> On Aug 10, 2012, at 3:29 PM, Scott Makeig <smakeig at gmail.com> wrote:
>>
>> MAX - Unfortunately, in general using two dipoles rather than one will
>> ~always improve the fit. Even if the source is a pure single dipole, a
>> second dipole can be used to correct for noise or errors in the forward
>> head model. This is less always the case for the constrained
>> spatially-symmetric dipole pairs allowed by dipfit(). However, we have not
>> thought of an optimal way to decide between using one or (occasionally) two
>> dipoles to fit e.g. maps of ICA brain sources. The goal would be to decide
>> whether the two-dipole version is fitting noise/forward model error vs
>> actual bilateral source generation...
>>
>> Scott Makeig
>>
>> On Thu, Aug 9, 2012 at 1:54 AM, Maximilien Chaumon <
>> maximilien.chaumon at gmail.com> wrote:
>>
>>> Hello all,
>>>
>>> When fitting dipoles to components, we are all sooner or later puzzled
>>> by the question whether to use one or two symmetrical dipoles.
>>>
>>> Would it be correct to put the problems in terms of a nested hypothesis
>>> testing?
>>>
>>> We are fitting a scalp map with one or two parameters and get a residual
>>> variance after the fit.
>>> Could we not use this residual variance as a measure of the SSE and
>>> compute a F statistic to decide whether to use the more complex (with two
>>> dipoles) or simpler (with one dipole) of two nested models?
>>> If yes, then how would we decide on the number of degrees of freedom?
>>> How many free parameters do we have in each case? x,y,z,and two
>>> orientations per dipole? how does the imposed symmetry affect that number?
>>> Could we really map residual variance to SSE? How many "observations" do we
>>> have in that case (see formula below)?
>>>
>>> I found this formula, for F:
>>> F = (SSEF-SSER)/ (kF-kR) / ((1-SSEF)/(N-kF-1))
>>> where
>>> SSE is sum of squared errors,
>>> k is numbers of parameters,
>>> N is number of observations (? what in our case?)
>>> F and R indices for full and reduced model respectively (in our case two
>>> and one dipole).
>>>
>>>
>>> Thanks a lot for any comment!
>>> Best,
>>> Max
>>>
>>> dipfit
>>>
>>>
>>>
>>>
>>>
>>> _______________________________________________
>>> Eeglablist page: http://sccn.ucsd.edu/eeglab/eeglabmail.html
>>> To unsubscribe, send an empty email to
>>> eeglablist-unsubscribe at sccn.ucsd.edu
>>> For digest mode, send an email with the subject "set digest mime" to
>>> eeglablist-request at sccn.ucsd.edu
>>>
>>
>>
>>
>> --
>> Scott Makeig, Research Scientist and Director, Swartz Center for
>> Computational Neuroscience, Institute for Neural Computation; Prof. of
>> Neurosciences (Adj.), University of California San Diego, La Jolla CA
>> 92093-0559, http://sccn.ucsd.edu/~scott
>>
>> _______________________________________________
>> Eeglablist page: http://sccn.ucsd.edu/eeglab/eeglabmail.html
>> To unsubscribe, send an empty email to
>> eeglablist-unsubscribe at sccn.ucsd.edu
>> For digest mode, send an email with the subject "set digest mime" to
>> eeglablist-request at sccn.ucsd.edu
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20120824/32860f16/attachment-0001.html>
More information about the fieldtrip
mailing list