fmri_perm_behav.m

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

function [brainlv,s,behavlv,brainscores,behavscores,lvcorrs, ...
	origpost,perm_result,behavdata,evt_list,datamatcorrs_lst]= ...
		fmri_perm_behav(datamat,behavdata,evt_list, ...
		behavdata_lst, newdata_lst, num_subj_lst, ...
		num_perm,num_cond,num_subj,posthoc)

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

    datamatcorrs_lst = {};

    stacked_datamatcorrs = [];
    stacked_datamat = datamat;
    stacked_behavdata = behavdata;

    perm_result = [];

    num_groups = length(newdata_lst);

    progress_hdl = rri_progress_ui('initialize');

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

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

       k = num_cond;
       n = num_subj_lst(i);

       datamat = newdata_lst{i};

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

       % compute correlation
       %
       datamatcorrs = rri_corr_maps(behavdata_lst{i}, datamat, n, k);

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

       stacked_datamatcorrs = [stacked_datamatcorrs; datamatcorrs];
       datamatcorrs_lst = [datamatcorrs_lst, {datamatcorrs}];

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

    end		% for

    % Singular Value Decomposition
    %
    [r c] = size(stacked_datamatcorrs);
    if r <= c
       [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
    %
    [brainscores, behavscores, lvcorrs] = ...
		rri_get_behavscores(stacked_datamat, stacked_behavdata, ...
		brainlv, behavlv, num_cond, num_subj_lst);

    rri_progress_ui(progress_hdl,'',1);

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

    for p = 1:num_perm
        reorder(:,p) = [randperm(size(stacked_datamat,1))'];
    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);

        % data_p = stacked_datamat(reorder(:,p),:);

        behav_p = stacked_behavdata(reorder(:,p),:);

        stacked_data = [];

        for g=1:num_groups

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

		% 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,:)));
                count = 0;
                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,:)));
                    count = count + 1;
                    if count > 100
                       msg = 'Please check your behavior data, and make ';
                       msg = [msg 'sure none of the columns are all the '];
                       msg = [msg 'same for each group'];
                       uiwait(msgbox(msg, 'Program can not proceed', 'modal'));
                       brainlv = [];
                       return;
                    end
                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

            stacked_data = [stacked_data; data];

        end		% for num_groups
	
        [r c] = size(stacked_data);
        if r <= c
           [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;					% fmri_perm_behav