fmri_perm_multiblock.m

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

function [brainlv,s,designlv,behavlv,brainscores,designscores,behavscores, ...
	lvcorrs,origpost,perm_result,behavdata,evt_list,behavdata_lst,datamatcorrs_lst, ...
	b_scores, behav_row_idx]= ...
		fmri_perm_multiblock(datamat,behavdata,evt_list, ...
		behavdata_lst, newdata_lst, num_subj_lst, ...
		num_perm,num_cond,num_subj,posthoc,bscan)

    % Init
    %
    brainlv = [];
    s = [];
    designlv = [];
    behavlv = [];
    brainscores = [];
    designscores = [];
    behavscores = [];
    lvcorrs = [];
    origpost = [];
    b_scores = [];
    datamatcorrs_lst = {};

    stacked_TBdatamatcorrs = [];
    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);

       k = num_cond;
       kk = length(bscan);

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

       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 task mean
       %
       Tdatamatcorrs = rri_task_mean(datamat,n)-ones(k,1)*mean(datamat);

       % compute correlation
       %
       Bdatamatcorrs = rri_corr_maps_notall(behavdata_lst{i}, datamat, n, bscan);

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

       %  stack task and behavior - keep un-normalize data that will be
       %  used to recover the normalized one
       %
       TBdatamatcorrs = [Tdatamatcorrs; Bdatamatcorrs];

       %  stack task and behavior - normalize to unit length to reduce
       %  scaling difference
       %
       datamatcorrs = [normalize(Tdatamatcorrs,2); normalize(Bdatamatcorrs,2)];

       stacked_TBdatamatcorrs = [stacked_TBdatamatcorrs; TBdatamatcorrs];
       stacked_datamatcorrs = [stacked_datamatcorrs; datamatcorrs];
       datamatcorrs_lst = [datamatcorrs_lst, {Bdatamatcorrs}];

         tmp_row_idx = reshape(1:size(stacked_behavdata,2)*k, [size(stacked_behavdata,2) k]);
         tmp_row_idx = tmp_row_idx(:,bscan);
         row_idx = [row_idx tmp_row_idx(:)+size(stacked_behavdata,2)*k*(i-1)];

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

    end		% for

    if ~isempty(posthoc)
       row_idx = row_idx(:);
       posthoc = posthoc(row_idx,:);
    end

    % Singular Value Decomposition
    %
    [r c] = size(stacked_datamatcorrs);
    if r <= c
       [brainlv, s, v] = svd(stacked_datamatcorrs',0);
    else
       [v, s, brainlv] = svd(stacked_datamatcorrs,0);
    end

    s = diag(s);
    original_v = v * diag(s);

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


   %  Since the 2 matrices that went into the SVD were unit normal, we should
   %  go backwards from the total Singular value Sum of Squares (SSQ)

      %  Calculate total SSQ
      %
      total_s = sum(stacked_TBdatamatcorrs(:).^2);

      %  Calculate distribution of normalized SSQ across LVs
      %
      per = s.^2 / sum(s.^2);

      %  Re-calculate singular value based on the distribution of SSQ
      %  across normalized LVs
      %
      org_s = sqrt(per * total_s);

      %  Re-scale v (block LV) with singular value
      %
      org_v = v * diag(org_s);


      %  Separate v into 2 parts: designlv and behavlv
      %
      designlv = [];
      behavlv = [];

      for g = 1:num_groups
         t = size(stacked_behavdata, 2);

         designlv = [designlv; v((g-1)*k+(g-1)*kk*t+1 : (g-1)*k+(g-1)*kk*t+k,:)];
         behavlv = [behavlv; v((g-1)*k+(g-1)*kk*t+k+1 : (g-1)*k+(g-1)*kk*t+k+kk*t,:)];
      end

      num_col = size(designlv, 2);

      % expand the num_subj for each row (cond)
      % did the samething as testvec
      %
      designscores = [];
      row_idx = [];
      last = 0;

      for g = 1:num_groups
         n = num_subj_lst(g);

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

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

         %  take this advantage (having g & n) to get row_idx
         %
         tmp = 1:n*k;
         tmp = reshape(tmp, [n k]);
         tmp = tmp(:, bscan);
         behavdata_lst{g} = behavdata_lst{g}(tmp(:),:);

         row_idx = [row_idx ; tmp(:) + last];
         last = last + n*k;
      end

    behav_row_idx = row_idx;
    behavdata = behavdata(row_idx,:);

    % calculate behav scores
    %
    b_scores = stacked_datamat * brainlv;
    [brainscores, behavscores, lvcorrs] = ...
		rri_get_behavscores(stacked_datamat(row_idx,:), ...
		stacked_behavdata(row_idx,:), ...
		brainlv, behavlv, kk, num_subj_lst);

%      lvcorrs = original_v;

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


    rri_progress_ui(progress_hdl,'',1);

    %  Begin permutation loop
    %
    sp = zeros(size(s));
    dp = zeros(size(v));

    Treorder = rri_perm_order(num_subj_lst, k, num_perm);
    rand('state',sum(100*clock));

    for p = 1:num_perm
        Breorder(:,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(Treorder(:,p),:);
        behav_p = stacked_behavdata(Breorder(:,p),:);

        stacked_TBdata = [];
        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;

                if num_groups == 1
                  Tdata = rri_task_mean(data_p,n)-ones(k,1)*mean(data_p);
                else
                  Tdata = rri_task_mean(data_p(1+span:n*k+span,:),n)-ones(k,1)*mean(data_p(1+span:n*k+span,:));
                end

                min1 = min(std(behav_p(1+span:n*k+span,:)));
                count = 0;
                while (min1 == 0)
                    Breorder(:,p) = [randperm(size(stacked_datamat,1))'];
                    behav_p = stacked_behavdata(Breorder(:,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
                   Bdata = rri_corr_maps_notall(behav_p, stacked_datamat, n, bscan);
                else
                   Bdata = rri_corr_maps_notall(behav_p(1+span:n*k+span,:), ...     
				stacked_datamat(1+span:n*k+span,:), n, bscan);
                end

            TBdata = [Tdata; Bdata];
            data = [normalize(Tdata,2); normalize(Bdata,2)];

            stacked_TBdata = [stacked_TBdata; TBdata];
            stacked_data = [stacked_data; data];

        end		% for num_groups
	
        [r c] = size(stacked_data);
        if r <= c
           [pbrainlv, sperm, pv] = svd(stacked_data',0);
        else
           [pv, sperm, pbrainlv] = svd(stacked_data,0);
        end

        rotatemat = rri_bootprocrust(v,pv);
        pv = pv * sperm * rotatemat;
        sperm = sqrt(sum(pv.^2));

        ptotal_s = sum(stacked_TBdata(:).^2);
        per = diag(sperm).^2 / sum(diag(sperm).^2);
        sperm = sqrt(per * ptotal_s);
        pv = normalize(pv) * diag(sperm);

        sp = sp + (sperm>=org_s);
        dp = dp + (abs(pv) >= abs(org_v));

        Bpv = [];

        for g = 1:num_groups
            t = size(stacked_behavdata, 2);
            Bpv = [Bpv; pv((g-1)*k+(g-1)*kk*t+k+1 : (g-1)*k+(g-1)*kk*t+k+kk*t,:)];
        end

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

    end		% for num_perm

    if num_perm ~= 0

        perm_result.sprob = sp ./ (num_perm + 1);
        perm_result.s_prob = perm_result.sprob;
        perm_result.vprob = dp ./ (num_perm + 1);
        perm_result.num_perm = num_perm;
        perm_result.Tpermsamp = Treorder;
        perm_result.Bpermsamp = Breorder;
        perm_result.sp = sp;
        % perm_result.dp = dp;

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

    end

    return;					% fmri_perm_multiblock