Hi everyone,<br><br>Just to follow-up on this if anyone wanted to take that approach.<br>I did use the following procedure to decide whether to fit one or two dipoles to a given component.<br>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))<br>
2) fit two dipoles, compute SSE for each component map: sse(2)<br>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<br>
4) F statistic : F = (sse(1)-sse(2))/(p2-p1)/(sse(2)/(n-p2))<br>5) if F > finv(.95,p2-p1,n-p2) decide to go for 2 dipoles.<br>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.<br>
<br>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.<br>
<br>One would probably better just screen all components by eye and spot those that look dipolar using thumbometric visual assessment.<br><br>Best,<br>Max<br><br>PS: the script I used:<br><br>dofit = 1;<br>dorework = 0;<br>
doplot = 1;<br><br>rootdir = '/DATA/';<br>cd(rootdir)<br>studyfile = fullfile(rootdir,'AP2.study');<br><br>if dofit || dorework<br> [STUDY ALLEEG] = pop_loadstudy(studyfile);<br> <br> for i_set = 1:numel(STUDY.datasetinfo)<br>
EEG = pop_loadset('filename',ALLEEG(i_set).filename,'filepath',ALLEEG(i_set).filepath);%pop_loadset(GRAB(i_set).name);<br> if isfield(EEG,'peakalphafrequency')<br> EEG.peakalphafreq = EEG.peakalphafrequency;<br>
EEG = rmfield(EEG,'peakalphafrequency');<br> end<br> <br> if ~dorework<br> 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] );<br>
EEGonedipole = pop_multifit(EEG, [1:size(EEG.icawinv,2)] ,'threshold',100,'rmout','on','dipoles',1);<br> EEGtwodipoles = pop_multifit(EEG, [1:size(EEG.icawinv,2)] ,'threshold',100,'rmout','on','dipoles',2);<br>
else<br> EEGonedipole = EEG;<br> EEGtwodipoles = EEG;<br> EEGonedipole.dipfit.model = EEG.dipfit.onedipole;<br> EEGtwodipoles.dipfit.model = EEG.dipfit.twodipoles;<br> end<br>
% H0: one dipole.<br> EEG.dipfit = EEGonedipole.dipfit;<br> % still we store them both. doesn't hurt.<br> EEG.dipfit.onedipole = EEGonedipole.dipfit.model;<br> EEG.dipfit.twodipoles = EEGtwodipoles.dipfit.model;<br>
% now for every component, if the rv is much better for the<br> % 2dipoles than for the 1dipole, we'll use the 2dipoles.<br> for i_dip = 1:numel(EEGonedipole.dipfit.model)<br> sse(1) = EEGonedipole.dipfit.model(i_dip).rv .* sum(EEGonedipole.icawinv(:,i_dip).^2);<br>
sse(2) = EEGtwodipoles.dipfit.model(i_dip).rv .* sum(EEGtwodipoles.icawinv(:,i_dip).^2);<br> p1 = 6;% six params for 1dipole model<br> p2 = 9;% 9 for 2 dipoles model (assume only symmetry is fixed).<br>
n = EEG.nbchan; % N is number of electrodes<br> <br> % F statistic<br> F = (sse(1)-sse(2))/(p2-p1)/(sse(2)/(n-p2));<br> EEG.dipfit.Fone2two(i_dip) = F;<br> % F decision<br>
Fthresh = finv(.999,p2-p1,n-p2);<br> if F > Fthresh % reject H0, take 2 dipoles model<br> % only if x ~= 0 and the two momenta are not just opposite.<br> mom = EEGtwodipoles.dipfit.model(i_dip).momxyz;<br>
pos = EEGtwodipoles.dipfit.model(i_dip).posxyz;<br> [t p r] = cart2sph(mom(:,1),mom(:,2),mom(:,3));<br> t(2) = t(2)+pi;%<br> d = sqrt(sum(diff(pos).^2));<br> % if the two dipoles are too close and have exactly<br>
% opposite angles<br> cond = d < 5 && (diff(abs(t)) < 0.01 || diff(abs(p)) < 0.01);<br> % or if they are exactly on the sagittal plane<br> cond = cond | abs(EEGtwodipoles.dipfit.model(i_dip).posxyz(1)) < 5;<br>
% or if the norm of one is much bigger than the other<br> ratio = norm(mom(1,:)) ./ norm(mom(2,:));<br> ratiothresh = 5;<br> cond = cond | (ratio > ratiothresh || ratio < 1/ratiothresh);<br>
% we will still keep the one dipole fit<br> if ~cond% so if not, we take the 2 dipole fit.<br> EEG.dipfit.model(i_dip) = EEGtwodipoles.dipfit.model(i_dip);<br> end<br>
end<br> end<br> <br> % reject poor overall fits.<br> EEG.dipfit.model = dipfit_reject(EEG.dipfit.model, 20/100);<br> % reject NaN rv<br> for i_m = 1:numel(EEG.dipfit.model)<br>
if isnan(EEG.dipfit.model(i_m).rv)<br> EEG.dipfit.model(i_m).posxyz = [];<br> EEG.dipfit.model(i_m).momxyz = [];<br> EEG.dipfit.model(i_m).rv = 1;<br> end<br>
end<br> <br> % removing dipoles outside the head<br> % ---------------------------------<br> rmdip = [];<br> for index = 1:numel(EEGonedipole.dipfit.model)<br> if ~isempty(EEG.dipfit.model(index).posxyz)<br>
if any(sqrt(sum(EEG.dipfit.model(index).posxyz.^2,2)) > 85)<br> rmdip = [ rmdip index];<br> EEG.dipfit.model(index).posxyz = [];<br> EEG.dipfit.model(index).momxyz = [];<br>
EEG.dipfit.model(index).rv = 1;<br> end;<br> end;<br> end;<br> if ~isempty(rmdip)<br> fprintf('%d out of cortex dipoles removed (usually artifacts)\n', length(rmdip));<br>
end;<br> <br> pop_saveset(EEG,'savemode','resave');<br> <br> end<br> <br>end<br>%%<br>if doplot<br> % [STUDY ALLEEG] = pop_loadstudy('AlphaPsyc2.study');<br> <br>
num = [];% number of dipoles that are fit for each component<br> idx = [];<br> topos = {[] [] []};<br> Fdiff = {[] [] []};<br> idx = {[] [] []};<br> dips = {};<br> for i_e = 1:numel(ALLEEG)<br> for i_comp = 1:numel(ALLEEG(i_e).dipfit.model)<br>
num(end+1) = size(ALLEEG(i_e).dipfit.model(i_comp).posxyz,1);<br> topos{num(end)+1}(end+1,:) = ALLEEG(i_e).icawinv(:,i_comp)';<br> dips{num(end)+1}(size(topos{num(end)+1},1)) = ALLEEG(i_e).dipfit.model(i_comp);<br>
Fdiff{num(end)+1}(end+1) = ALLEEG(i_e).dipfit.Fone2two(i_comp);<br> idx{num(end)+1}(end+1,:) = [i_e i_comp];<br> end<br> end<br> for i_plot = 2:3<br> figure(6548+i_plot);paste_figpos manytopos<br>
ntopos = size(topos{i_plot},1);<br> n = ceil(sqrt(ntopos));<br> m = n;<br> for i_topo = 1:ntopos<br> subplot(n,m,i_topo)<br> cla<br> topoplot(topos{i_plot}(i_topo,:),ALLEEG(1).chanlocs,'conv','on','electrodes','off'); <br>
title({num2str(idx{i_plot}(i_topo,:)) num2str(Fdiff{i_plot}(i_topo))});<br> drawnow<br> end<br> end<br> <br>end<br><br><br><br><br><div class="gmail_quote">2012/8/13 Maximilien Chaumon <span dir="ltr"><<a href="mailto:maximilien.chaumon@gmail.com" target="_blank">maximilien.chaumon@gmail.com</a>></span><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hello eeglab & Fieldtrip,<br><br>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.<br>
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.<br>
<br>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.<br>Christian then said it'd be worth a shot, and I agree, so here I am again with two questions, or two confirmations:<br>
<br>1) <b>How many parameters are estimated in ft_dipolefitting.m ?</b> 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.<br>
2) <b>Can I assimilate the relative residual variance to a SSE?</b> the function rv.m does this: rv = sum((d1-d2).^2) ./ sum(d1.^2);<br>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?<br>
<br>Thanks a lot!<br>Max<div class="HOEnZb"><div class="h5"><br><br><div class="gmail_quote">
2012/8/11 Christian Kothe <span dir="ltr"><<a href="mailto:christiankothe@gmail.com" target="_blank">christiankothe@gmail.com</a>></span><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div bgcolor="#FFFFFF"><div>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).</div>
<div><br></div><div>The number of parameters for a 2-dipole model seems to be 3 (xyz) +<span> 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.</span></div>
<span><font color="#888888"><div><br>Christian</div></font></span><div><div><div><br>On Aug 10, 2012, at 3:29 PM, Scott Makeig <<a href="mailto:smakeig@gmail.com" target="_blank">smakeig@gmail.com</a>> wrote:<br>
<br></div><div></div><blockquote type="cite"><div>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...<div>
<br></div><div>Scott Makeig<br><br><div class="gmail_quote">On Thu, Aug 9, 2012 at 1:54 AM, Maximilien Chaumon <span dir="ltr"><<a href="mailto:maximilien.chaumon@gmail.com" target="_blank">maximilien.chaumon@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hello all,<br><br>When fitting dipoles to components, we are all sooner or later puzzled by the question whether to use one or two symmetrical dipoles.<br>
<br>Would it be correct to put the problems in terms of a nested hypothesis testing?<br>
<br>We are fitting a scalp map with one or two parameters and get a residual variance after the fit.<br>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?<br>
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)?<br>
<br>I found this formula, for F:<br>F = (SSEF-SSER)/ (kF-kR) / ((1-SSEF)/(N-kF-1))<br>where<br>SSE is sum of squared errors,<br>k is numbers of parameters,<br>N is number of observations (? what in our case?)<br>F and R indices for full and reduced model respectively (in our case two and one dipole).<br>
<br><br>Thanks a lot for any comment!<br>Best,<br>Max<br><br>dipfit<br><br><br><br><br>
<br>_______________________________________________<br>
Eeglablist page: <a href="http://sccn.ucsd.edu/eeglab/eeglabmail.html" target="_blank">http://sccn.ucsd.edu/eeglab/eeglabmail.html</a><br>
To unsubscribe, send an empty email to <a href="mailto:eeglablist-unsubscribe@sccn.ucsd.edu" target="_blank">eeglablist-unsubscribe@sccn.ucsd.edu</a><br>
For digest mode, send an email with the subject "set digest mime" to <a href="mailto:eeglablist-request@sccn.ucsd.edu" target="_blank">eeglablist-request@sccn.ucsd.edu</a><br></blockquote></div><br><br clear="all">
<div><br></div>
-- <br>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, <a href="http://sccn.ucsd.edu/%7Escott" target="_blank">http://sccn.ucsd.edu/~scott</a><br>
</div>
</div></blockquote><blockquote type="cite"><div><span>_______________________________________________</span><br><span>Eeglablist page: <a href="http://sccn.ucsd.edu/eeglab/eeglabmail.html" target="_blank">http://sccn.ucsd.edu/eeglab/eeglabmail.html</a></span><br>
<span>To unsubscribe, send an empty email to <a href="mailto:eeglablist-unsubscribe@sccn.ucsd.edu" target="_blank">eeglablist-unsubscribe@sccn.ucsd.edu</a></span><br><span>For digest mode, send an email with the subject "set digest mime" to <a href="mailto:eeglablist-request@sccn.ucsd.edu" target="_blank">eeglablist-request@sccn.ucsd.edu</a></span></div>
</blockquote></div></div></div></blockquote></div><br>
</div></div></blockquote></div><br>