pop_comperp.m

含有多种ICA算法的eeglab工具箱

srate  = ALLEEG(datadd(1)).srate;
xmin   = ALLEEG(datadd(1)).xmin;
xmax   = ALLEEG(datadd(1)).xmax;
nbchan = ALLEEG(datadd(1)).nbchan;
chanlocs = ALLEEG(datadd(1)).chanlocs;
for index = union(datadd, datsub)
    if ALLEEG(index).pnts ~= pnts,     error(['Dataset '  int2str(index) ' does not have the same number of points as others']); end;
    if ALLEEG(index).xmin ~= xmin,     error(['Dataset '  int2str(index) ' does not have the same xmin as others']); end;
    if ALLEEG(index).xmax ~= xmax,     error(['Dataset '  int2str(index) ' does not have the same xmax as others']); end;
    if ALLEEG(index).nbchan ~= nbchan, error(['Dataset '  int2str(index) ' does not have the same number of channels as others']); end;
end;
    
% compute ERPs for add
% --------------------
for index = 1:length(datadd)
    TMPEEG = eeg_checkset(ALLEEG(datadd(index)));
    if flag == 1, erp1ind(:,:,index)  = mean(TMPEEG.data,3);
    elseif isempty(TMPEEG.icaact)
        tmpica              =  reshape((TMPEEG.icaweights*TMPEEG.icasphere)*TMPEEG.data(:,:), ...
                                       size(TMPEEG.icaweights,1), size(TMPEEG.data,2), size(TMPEEG.data,3));
        erp1ind(:,:,index)  = mean(tmpica,3);
    else          
        tmpica              =  reshape((TMPEEG.icaweights*TMPEEG.icasphere)*TMPEEG.data(:,:), ...
                                       size(TMPEEG.icaweights,1), size(TMPEEG.data,2), size(TMPEEG.data,3));
        erp1ind(:,:,index)  = mean(tmpica,3);
        erp1ind(:,:,index)  = mean(TMPEEG.icaact,3);
    end;
    addnames{index} = [ '#' int2str(datadd(index)) ' ' TMPEEG.setname ' (n=' int2str(TMPEEG.trials) ')' ];
    clear TMPEEG;
end;

% optional: subtract
% ------------------
colors = {}; % color aspect for curves
allcolors = { 'b' 'r' 'g' 'c' 'm' 'y' [0 0.5 0] [0.5 0 0] [0 0 0.5] [0.5 0.5 0] [0 0.5 0.5] [0.5 0 0.5] [0.5 0.5 0.5] };
allcolors = { allcolors{:} allcolors{:} allcolors{:} allcolors{:} allcolors{:} allcolors{:} };
allcolors = { allcolors{:} allcolors{:} allcolors{:} allcolors{:} allcolors{:} allcolors{:} };
if length(datsub) > 0 % dataset to subtract

    % compute ERPs for sub
    % --------------------
    for index = 1:length(datsub)
        TMPEEG = eeg_checkset(ALLEEG(datsub(index)));
        if flag == 1, erp2ind(:,:,index)  = mean(TMPEEG.data,3);
        elseif isempty(TMPEEG.icaact)
        else          tmpica              =  reshape((TMPEEG.icaweights*TMPEEG.icasphere)*TMPEEG.data(:,:), ...
                                                     size(TMPEEG.icaweights,1), size(TMPEEG.data,2), size(TMPEEG.data,3));
                      erp2ind(:,:,index)  = mean(tmpica,3);
                      erp2ind(:,:,index)  = mean(TMPEEG.icaact,3);
        end;
        subnames{index} = [ '#' int2str(datsub(index)) ' ' TMPEEG.setname '(n=' int2str(TMPEEG.trials) ')' ];
        clear TMPEEG
    end;
    
    l1 = size(erp1ind,3);
    l2 = size(erp2ind,3);
    allcolors1 = allcolors(3:l1+2);
    allcolors2 = allcolors(l1+3:l1+l2+3);
    allcolors3 = allcolors(l1+l2+3:end);
    [erps1, erpstd1, colors1, colstd1, legend1] = preparedata( erp1ind        , g.addavg , g.addstd , g.addall , g.mode, 'Add ' , addnames, 'b', allcolors1 );
    [erps2, erpstd2, colors2, colstd2, legend2] = preparedata( erp2ind        , g.subavg , g.substd , g.suball , g.mode, 'Sub ' , subnames, 'r', allcolors2 );
    [erps3, erpstd3, colors3, colstd3, legend3] = preparedata( erp1ind-erp2ind, g.diffavg, g.diffstd, g.diffall, g.mode, 'Diff ', ...
                                                      { addnames subnames }, 'k', allcolors3 );
    
    % handle special case of std
    % --------------------------
    erptoplot  = [ erps1 erps2 erps3 erpstd1 erpstd2 erpstd3 ];
    colors     = { colors1{:} colors2{:} colors3{:} colstd1{:} colstd2{:} colstd3{:}};
    legend     = { legend1{:} legend2{:} legend3{:} };
    
    % highlight significant regions
    % -----------------------------
    if ~isempty(g.alpha)
        pvalues = pttest(erp1ind, erp2ind, 3);
        regions = p2regions(pvalues, g.alpha, [xmin xmax]*1000);
    else 
        pvalues= [];
    end;
    
else
    [erptoplot, erpstd, colors, colstd, legend] = preparedata( erp1ind, g.addavg, g.addstd, g.addall, g.mode, '', addnames, 'k', allcolors);
    erptoplot = [ erptoplot erpstd ];
    colors    = { colors{:} colstd{:} };
    
    % highlight significant regions
    % -----------------------------
    if ~isempty(g.alpha)
        pvalues = ttest(erp1ind, 0, 3);
        regions = p2regions(pvalues, g.alpha, [xmin xmax]*1000);
    else 
        pvalues= [];
    end;
end;
    
% lowpass data
% ------------
if ~isempty(g.lowpass)
    if exist('filtfilt') == 2
        erptoplot = eegfilt(erptoplot, srate, 0, g.lowpass);
    else
        erptoplot = eegfiltfft(erptoplot, srate, 0, g.lowpass);
    end;
end;
if strcmpi(g.geom, 'array') | flag == 0, chanlocs = []; end;

% select time range
% -----------------
if ~isempty(g.tlim)
    pointrange = round(eeg_lat2point(g.tlim/1000, [1 1], srate, [xmin xmax]));
    g.tlim     = eeg_point2lat(pointrange, [1 1], srate, [xmin xmax]);
    erptoplot  = reshape(erptoplot, size(erptoplot,1), pnts, size(erptoplot,2)/pnts);
    erptoplot  = erptoplot(:,[pointrange(1):pointrange(2)],:);
    pnts       = size(erptoplot,2);
    erptoplot  = reshape(erptoplot, size(erptoplot,1), pnts*size(erptoplot,3));    
    xmin = g.tlim(1);
    xmax = g.tlim(2);
end;

% plot data
% ---------
plottopo( erptoplot, 'chanlocs', chanlocs, 'frames', pnts, ...
          'limits', [xmin xmax 0 0]*1000, 'title', g.title, 'colors', colors, ...
          'chans', g.chans, 'legend', legend, 'regions', regions, 'ylim', g.ylim, g.tplotopt{:});

% outputs
% -------
times  = linspace(xmin, xmax, pnts);
erp1   = mean(erp1ind,3);
if length(datsub) > 0 % dataset to subtract
    erp2   = mean(erp2ind,3);
    erpsub = erp1-erp2;
else
    erp2 = [];
    erpsub = [];
end;

if nargin < 3 & nargout == 1
    erp1 = sprintf('pop_comperp( %s, %d, %s);', inputname(1), ...
                  flag, vararg2str({ datadd datsub options{:} }) );
end;
return;

% convert significance values to alpha
% ------------------------------------
function regions = p2regions( pvalues, alpha, limits);
    
    for index = 1:size(pvalues,1)
        signif = diff([1 pvalues(index,:) 1] < alpha);
        pos    = find([signif] > 0);
        pos    = pos/length(pvalues)*(limits(2) - limits(1))+limits(1);
        neg    = find([signif(2:end)] < 0);
        neg    = neg/length(pvalues)*(limits(2) - limits(1))+limits(1);
        if length(pos) ~= length(neg), signif, pos, neg, error('Region error'); end;
        regions{index} = [neg;pos];
    end;
    
% process data
% ------------
function [erptoplot, erpstd, colors, colstd, legend] = preparedata( erpind, plotavg, plotstd, plotall, mode, tag, dataset, coloravg, allcolors);

    colors    = {};
    legend    = {};
    erptoplot = [];
    erpstd    = [];
    colstd    = {};

    % plot individual differences
    % ---------------------------
    if strcmpi(plotall, 'on')
        erptoplot = [ erptoplot erpind(:,:) ];
        for index=1:size(erpind,3)
            if iscell(dataset)
                if strcmpi(tag, 'Diff ')
                    legend = { legend{:} [ dataset{1}{index} ' - ' dataset{2}{index} ] };
                else
                    legend = { legend{:} dataset{index} };
                end;
            else
                legend = { legend{:} [ 'Dataset ' int2str(dataset(index)) ] };
            end;
            colors = { colors{:}  allcolors{index} };
        end;
    end;
    
    % plot average
    % ------------
    if strcmpi( plotavg, 'on')
        if strcmpi(mode, 'ave')
             granderp    = mean(erpind,3);
             legend      = { legend{:} [ tag 'Average' ] };
        else granderp    = sqrt(mean(erpind.^2,3));
             legend      = { legend{:} [ tag 'RMS' ] };
        end;
        colors    = { colors{:}  {coloravg 'linewidth' 2 }};
        erptoplot = [ erptoplot granderp];
    end;

    % plot standard deviation
    % -----------------------
    if strcmpi(plotstd, 'on')
        if strcmpi(plotavg, 'on')
            std1      = std(erpind, [], 3);
            erptoplot = [ erptoplot granderp+std1 ];
            erpstd    = granderp-std1;
            legend    = { legend{:} [ tag 'Standard dev.' ] };
            colors    = { colors{:} { coloravg 'linestyle' ':' } };
            colstd    = { { coloravg 'linestyle' ':' } };
        else 
            disp('Warning: cannot show standard deviation without showing average');
        end;
    end;

% ------------------------------------------------------------------
    
function [p, t, df] = pttest(d1, d2, dim)
%PTTEST Student's paired t-test.
%       PTTEST(X1, X2) gives the probability that Student's t
%       calculated on paired data X1 and X2 is higher than
%       observed, i.e. the "significance" level. This is used
%       to test whether two paired samples have significantly
%       different means.
%       [P, T] = PTTEST(X1, X2) gives this probability P and the
%       value of Student's t in T. The smaller P is, the more
%       significant the difference between the means.
%       E.g. if P = 0.05 or 0.01, it is very likely that the
%       two sets are sampled from distributions with different
%       means.
%
%       This works for PAIRED SAMPLES, i.e. when elements of X1
%       and X2 correspond one-on-one somehow.
%       E.g. residuals of two models on the same data.
%       Ref: Press et al. 1992. Numerical recipes in C. 14.2, Cambridge.

if size(d1,dim) ~= size(d2, dim)
   error('PTTEST: paired samples must have the same number of elements !')
end
if size(d1,dim) == 1
    close; error('Cannot compute paired t-test for a single ERP difference')
end; 
a1 = mean(d1, dim);
a2 = mean(d2, dim);
v1 = std(d1, [], dim).^2;
v2 = std(d2, [], dim).^2;
n1 = size(d1,dim);
df = n1 - 1;
disp(['Computing t-values, df:' int2str(df) ]);

d1 = d1-repmat(a1, [ones(1,dim-1) size(d1,3)]);
d2 = d2-repmat(a2, [ones(1,dim-1) size(d2,3)]);
%cab = (x1 - a1)' * (x2 - a2) / (n1 - 1);
cab = mean(d1.*d2,3)/(n1-1);
% use abs to avoid numerical errors for very similar data
% for which v1+v2-2cab may be close to 0.
t = (a1 - a2) ./ sqrt(abs(v1 + v2 - 2 * cab) / n1) ;
p = betainc( df ./ (df + t.*t), df/2, 0.5) ;

% ------------------------------------------------------------------

function [p, t] = ttest(d1, d2, dim)
%TTEST Student's t-test for equal variances.
%       TTEST(X1, X2) gives the probability that Student's t
%       calculated on data X1 and X2, sampled from distributions
%       with the same variance, is higher than observed, i.e.
%       the "significance" level. This is used to test whether
%       two sample have significantly different means.
%       [P, T] = TTEST(X1, X2) gives this probability P and the
%       value of Student's t in T. The smaller P is, the more
%       significant the difference between the means.
%       E.g. if P = 0.05 or 0.01, it is very likely that the
%       two sets are sampled from distributions with different
%       means.
%
%       This works if the samples are drawn from distributions with
%       the SAME VARIANCE. Otherwise, use UTTEST.
%
%See also: UTTEST, PTTEST.
if size(d1,dim) == 1
    close; error('Cannot compute t-test for a single ERP')
end; 
a1 = mean(d1, dim);
v1 = std(d1, [], dim).^2;
n1 = size(d1,dim);
if length(d2) == 1 & d2 == 0
    a2 = 0;
    n2 = n1;
    df = n1 + n2 - 2;
    pvar = (2*(n1 - 1) * v1) / df ;
else 
    a2 = mean(d2, dim);
    v2 = std(d2, [], dim).^2;
    n2 = size(d2,dim);
    df = n1 + n2 - 2;
    pvar = ((n1 - 1) * v1 + (n2 - 1) * v2) / df ;     
end;
disp(['Computing t-values, df:' int2str(df) ]);

t = (a1 - a2) ./ sqrt( pvar * (1/n1 + 1/n2)) ;
p = betainc( df ./ (df + t.*t), df/2, 0.5) ;