rri_analysis_perm.m

绝对经典,老外制作的功能强大的matlab实现PLS_TOOBOX

%RRI_ANALYSIS_PERM Apply Behavioral or Task PLS test and permutation test
%	on RRI scan.
%
%   Usage: [brainlv,s,behavlv,brainscores,behavscores,lvcorrs, ...
%                origpost, perm_result] = ...
%                rri_analysis_perm(ismean, ishelmert, iscontrast, isbehav, ...
%                newdata_lst,num_cond_lst,num_subj_lst, ...
%                behavdata_lst, helmertdata_lst, contrastdata_lst, ...
%                num_perm, posthoc);
%
%   See also PLS_PERM_TEST, PLS_DEVIATION_PERM_TEST, BEHAVPLS_PERM, TASKPLS_PERM
%

%   I (ismean) - 1 if using grand mean deviation;
%   I (ishelmert) - 1 if using helmert matrix;
%   I (iscontrast) - 1 if using contrast data;
%   I (isbehav) - 1 if using behavior data;
%   I (newdata_lst) - A group list of datamat files;
%   I (num_cond_lst) - A group list of condition numbers;
%   I (num_subj_lst) - A group list of subject numbers;
%   I (behavdata_lst) - A group list of behav data with selected columns;
%   I (helmertdata_lst) - A group list of helmert matrix with selected columns;
%   I (contrastdata_lst) - A group list of contrast data with selected columns;
%			contrastdata_lst has been orthchecked in
%			erp_get_common function
%   I (num_perm) - Number of Permutation;
%   I (posthoc) - posthoc data;
%
%   O (brainlv) - Left singular value vector. It is LV for the brain dimension.
%   O (s) - Singular values vector.
%   O (behavlv) - Right singular vector. It is LV for the behavior OR contrast
%		dimension for each scan.
%   O (brainscores) - The brain score.
%   O (behavscores) - If use GRAND_MEAN method, behavscores = behavlv, then
%		expand each condition for all the subjects. If not, it is
%		calculated for each set of condition, then stack together.
%   O (lvcorrs) - Correlates brain scores with behavior data for multiple scans,
%		return [] if for Task PLS or if use GRAND_MEAN method.
%   O (origpost) - posthoc result.
%   O (perm_result) - A Structure array containing the permutation result data.
%
%   Created on 03-OCT-2002 by Jimmy Shen for PLS test
%   Modified on 24-OCT-2002 by Jimmy Shen to add permutation test
%   Modified Jan 10,03 to add contrast & helmert
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [brainlv,s,behavlv,brainscores,behavscores,lvcorrs, ...
                origpost, perm_result, datamatcorrs_lst] = ...
                rri_analysis_perm(ismean, ishelmert, iscontrast, isbehav, ...
                newdata_lst,num_cond_lst,num_subj_lst, ...
                behavdata_lst, helmertdata_lst, contrastdata_lst, ...
                num_perm, posthoc);

    % Init
    %
    brainlv = [];
    s = [];
    behavlv = [];
    brainscores = [];
    behavscores = [];
    lvcorrs = [];

    datamatcorrs_lst = {};
    single_cond_lst = {};

    stacked_behavdata = [];
    stacked_helmertdata = [];
    stacked_contrastdata = [];
    stacked_datamatcorrs = [];
    stacked_datamat = [];
    stacked_mask = [];

    perm_result = [];

    num_groups = length(newdata_lst);

    progress_hdl = rri_progress_ui('initialize');

    msg = 'Working on PLS ...';
    rri_progress_ui(progress_hdl, '', msg);

    if ~isbehav & num_cond_lst(1)==1
       tmp = [];

       for g = 1:num_groups
          tmp = [tmp; newdata_lst{g}];
       end

       single_cond_lst = {tmp};
    end

    % loop accross the groups, and
    % calculate datamatcorrs for each group
    %
    for i = 1:num_groups

        k = num_cond_lst(i);
        n = num_subj_lst(i);

        if isempty(single_cond_lst)
           datamat = newdata_lst{i};
        elseif i==1
           datamat = single_cond_lst{1};
        end

        rri_progress_ui(progress_hdl,'',2/10+5/10*(i-1)/(num_groups)+1/(10*num_groups));

        % compute correlation or covariance
        %
        if ismean
           if isempty(single_cond_lst)
              datamatcorrs = rri_task_mean(datamat,n)-ones(k,1)*mean(datamat);
           elseif i==1
              datamatcorrs = rri_task_mean1(datamat,num_subj_lst)-ones(num_groups,1)*mean(datamat);
           end
        elseif ishelmert
            datamatcorrs = rri_xcovy(helmertdata_lst{i}, datamat);
        elseif iscontrast
            datamatcorrs = rri_xcovy(contrastdata_lst{i}, datamat);
        elseif isbehav
            datamatcorrs = rri_corr_maps(behavdata_lst{i}, datamat, n, k);
            datamatcorrs_lst = [datamatcorrs_lst, {datamatcorrs}];
        end

        rri_progress_ui(progress_hdl,'',2/10+5/10*(i-1)/(num_groups)+3/(10*num_groups));

	% if more than one group, stack data together
	%
        if ishelmert
            stacked_helmertdata = [stacked_helmertdata; helmertdata_lst{i}];
        elseif iscontrast
            stacked_contrastdata=[stacked_contrastdata;contrastdata_lst{i}];
        elseif isbehav
            stacked_behavdata = [stacked_behavdata; behavdata_lst{i}];
        end

        if isempty(single_cond_lst) | i==1
           stacked_datamat = [stacked_datamat; datamat];
           stacked_datamatcorrs = [stacked_datamatcorrs; datamatcorrs];
        end

        rri_progress_ui(progress_hdl,'',2/10+5/10*(i-1)/(num_groups)+5/(10*num_groups));

    end		% for

    % actually, all the groups must have the same condition number
    %
    num_cond = num_cond_lst(1);

    % Singular Value Decomposition
    %
    [r c] = size(stacked_datamatcorrs);
    if r <= c
        % transpose datamatcorrs to ensure SVD operation will be
        % on smallest of RxC dimension
        %
        [brainlv,s,behavlv] = svd(stacked_datamatcorrs',0);
    else
        [behavlv,s,brainlv] = svd(stacked_datamatcorrs,0);
    end

    s = diag(s);

    if ~isempty(posthoc)
        origpost = rri_xcor(posthoc,behavlv);
        porigpost = zeros(size(origpost));
    else
        origpost = [];
    end

    rri_progress_ui(progress_hdl,'',9/10);

    % calculate behav scores
    %
    if ismean

        brainscores = stacked_datamat * brainlv;

        % Here, behavlv is actually designlv
        % according to taskpls.m: fscores=design=testvec(designlv)
	% so, designscores = designlv

        num_col = size(behavlv, 2);

        % expand the num_subj for each row (cond)
        % did the samething as testvec
        %
        for i = 1:num_groups

            k = num_cond_lst(i);
            n = num_subj_lst(i);

            tmp = reshape(behavlv((i-1)*k+1:(i-1)*k+k,:), ...
		[1, num_col*k]);
            tmp = repmat(tmp, [n, 1]);		% expand to num_subj
            tmp = reshape(tmp, [n*k, num_col]);

            behavscores = [behavscores; tmp];	% stack by groups

        end

    elseif ishelmert | iscontrast
        if ishelmert
            behav = stacked_helmertdata;
        elseif iscontrast
            behav = stacked_contrastdata;
        end

        [brainscores, behavscores] = ...
		rri_get_designscores(stacked_datamat, behav, ...
		brainlv, behavlv, num_cond_lst(1), num_subj_lst);
    elseif isbehav
        behav = stacked_behavdata;

        [brainscores, behavscores, lvcorrs] = ...
		rri_get_behavscores(stacked_datamat, behav, ...
		brainlv, behavlv, num_cond_lst(1), num_subj_lst);
    end		% if

    rri_progress_ui(progress_hdl,'',1);

    %  Begin permutation loop
    %
    sp = zeros(size(s));
    dp = zeros(size(behavlv));
    rand('state',sum(100*clock));

    if isbehav
        for p = 1:num_perm
            reorder(:,p) = [randperm(size(stacked_datamat,1))'];
        end
    else
%        reorder = rri_mkperm_order(num_cond_lst(1), num_subj_lst, num_perm);
        reorder = rri_perm_order(num_subj_lst,num_cond_lst(1),num_perm);
    end

    for p = 1:num_perm

        msg = ['Working on Permutation:  ',num2str(p),' out of ',num2str(num_perm)];
        rri_progress_ui(progress_hdl, '', msg);
        rri_progress_ui(progress_hdl,'',p/num_perm);

        if ~isbehav
            data_p = stacked_datamat(reorder(:,p),:);
        end

        if ishelmert
            behav_p = stacked_helmertdata(reorder(:,p),:);
        elseif iscontrast
            behav_p = stacked_contrastdata(reorder(:,p),:);
        elseif isbehav
            behav_p = stacked_behavdata(reorder(:,p),:);
        end
        stacked_data = [];

        for g=1:num_groups

            k = num_cond_lst(g);
            n = num_subj_lst(g);
            span = sum(num_subj_lst(1:g-1)) * num_cond;

            if ismean
               if isempty(single_cond_lst)
                  if num_groups == 1
%                     meanmat = rri_task_mean(data_p, n);
 %                    data = meanmat - (ones(k,1)*mean(meanmat));
                     data = rri_task_mean(data_p,n)-ones(k,1)*mean(data_p);
                  else
%                     meanmat = rri_task_mean(data_p(1+span:n*k+span,:), n);
 %                    data = meanmat - (ones(k,1)*mean(meanmat));
                     data = rri_task_mean(data_p(1+span:n*k+span,:),n)-ones(k,1)*mean(data_p(1+span:n*k+span,:));
                  end
               elseif g==1
                  data = rri_task_mean1(data_p,num_subj_lst)-ones(num_groups,1)*mean(data_p);
               end
            elseif ishelmert
                if num_groups == 1
                   data = rri_xcovy(behav_p, data_p);
                else
                   data = rri_xcovy(behav_p(1+span:n*k+span,:), ...
				data_p(1+span:n*k+span,:));
                end
            elseif iscontrast
                if num_groups == 1
                   data = rri_xcovy(behav_p, data_p);
                else
                   data = rri_xcovy(behav_p(1+span:n*k+span,:), ...
				data_p(1+span:n*k+span,:));
                end
            elseif isbehav

		% Check for upcoming NaN and re-sample if necessary.
		% this only happened on behavior analysis, because the
		% 'xcor' inside of 'rri_corr_maps' contains a 'stdev', which
		% is a divident. If it is 0, it will cause divided by 0
		% problem.
		% since this happend very rarely, so the speed will not
		% be affected that much.
		%
                min1 = min(std(behav_p(1+span:n*k+span,:)));
                while (min1 == 0)
                    reorder(:,p) = [randperm(size(stacked_datamat,1))'];
                    behav_p = stacked_behavdata(reorder(:,p),:);
                    min1 = min(std(behav_p(1+span:n*k+span,:)));
                end

		% Notice here that stacked_datamat is used, instead of
		% boot_p. This is only for behavpls_perm.
		%
                if num_groups == 1
                   data = rri_corr_maps(behav_p, stacked_datamat, n, k);
                else
                   data = rri_corr_maps(behav_p(1+span:n*k+span,:), ...
				stacked_datamat(1+span:n*k+span,:), n, k);
                end
            end

            if isempty(single_cond_lst) | g==1
               stacked_data = [stacked_data; data];
            end

        end		% for num_groups

        % Singular Value Decomposition
        %
        [r c] = size(stacked_data);
        if r <= c
            % transpose datamatcorrs to ensure SVD operation will be
            % on smallest of RxC dimension
            %
            [pbrainlv, sperm, pbehavlv] = svd(stacked_data',0);
        else
            [pbehavlv, sperm, pbrainlv] = svd(stacked_data,0);
        end

        rotatemat = rri_bootprocrust(behavlv,pbehavlv);
        pbehavlv = pbehavlv * sperm * rotatemat;
        sperm = sqrt(sum(pbehavlv.^2));
        sp = sp + (sperm'>=s);
        dp = dp + (abs(pbehavlv) >= abs(behavlv));

        if ~isempty(posthoc)
            tmp = rri_xcor(posthoc, pbehavlv);
            porigpost = porigpost + (abs(tmp) >= abs(origpost));
        end

    end		% for num_perm

    if num_perm ~= 0

        perm_result.sprob = sp ./ num_perm;
        % perm_result.dprob = dp ./ num_perm;
        perm_result.num_perm = num_perm;
        perm_result.permsamp = reorder;
        perm_result.sp = sp;
        % perm_result.dp = dp;

        if ~isempty(posthoc)
            perm_result.posthoc_prob = porigpost / num_perm;
        end

    end

    return;					% rri_analysis_perm


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   load_contrast_file: turn the contrast file into a contrast matrix,
%		which is used for behavdata or design data.
%
%   I (contrast_file) - filespec of the contrast file
%   I (num_subj) - number of subjects in each condition
%   O (contrasts) - contrast matrix
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function contrasts = load_contrast_file(contrast_file, num_subj)

    load(contrast_file);

    num_contrasts = length(pls_contrasts);
    num_conditions = length(pls_contrasts(1).value);

    contrasts = zeros(num_conditions,num_contrasts);
    for i=1:num_contrasts
        contrasts(:,i) = pls_contrasts(i).value';
    end

    % expand each row of condition for all the subjects
    %
    tmp = contrasts(:)';
    tmp = repmat(tmp,num_subj,1);
    contrasts = reshape(tmp, num_subj*num_conditions, num_contrasts);

    return;						% load_contrast_file