pet_analysis_boot.m

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

%PET_ANALYSIS_BOOT Apply Behavioral or Task PLS test and bootstrap test
%	on PET scan.
%
%   Usage: [brainlv,s,behavlv,brainscores,behavscores,lvcorrs,boot_result, ...
%		datamatcorrs_lst] = ...
%		pet_analysis_boot(behavdata_lst,newdata_lst,num_cond_lst,...
%		num_subj_lst,num_boot,isbehav,Clim)
%
%   See also PLS_BOOT_TEST, PLS_DEVIATION_BOOT_TEST, BEHAVPLS_BOOT, TASKPLS_BOOT
%

%   Called by pet_analysis
%
%   I (behavdata_lst) - A group list of behav data with selected columns;
%   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 (num_boot) - Number of Bootstrap;
%   I (isbehav) - 1 if applying Behavioral PLS test; 0 if applying Task PLS test.
%   I (Clim) - upper limit of confidence interval estimated.
%   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 (boot_result) - A Structure array containing the bootstrap result data.
%
%   Created on 03-OCT-2002 by Jimmy Shen for PLS test
%   Modified on 27-OCT-2002 by Jimmy Shen to add bootstrap test
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [brainlv,s,behavlv,brainscores,behavscores,lvcorrs, ...
		boot_result, datamatcorrs_lst] = ...
		pet_analysis_boot(behavdata_lst,newdata_lst,num_cond_lst,...
		num_subj_lst,num_boot,isbehav,Clim, ...
		min_subj_per_group,is_boot_samples,boot_samples,new_num_boot)

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

    datamatcorrs_lst = {};
    single_cond_lst = {};

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

    num_groups = length(newdata_lst);

    % keeps track of number of times a new bootstrap had to be generated
    %
    countnewtotal=0;
    num_LowVariability_behav_boots = [];
    badbeh = [];


%    boot_progress = rri_progress_ui('initialize');

    if isbehav
%       [reorder, new_num_boot] = rri_boot_order(num_subj_lst, num_cond_lst(1), num_boot, 0, boot_progress, ...
       [reorder, new_num_boot] = rri_boot_order(num_subj_lst, num_cond_lst(1), num_boot, 0, ...
		min_subj_per_group,is_boot_samples,boot_samples,new_num_boot);
    else
%       [reorder, new_num_boot] = rri_boot_order(num_subj_lst, num_cond_lst(1), num_boot, 1, boot_progress, ...
       [reorder, new_num_boot] = rri_boot_order(num_subj_lst, num_cond_lst(1), num_boot, 1, ...
		min_subj_per_group,is_boot_samples,boot_samples,new_num_boot);
    end

    if isempty(reorder)
       return;
    end;


    progress_hdl = rri_progress_ui('initialize');

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


%    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 & covariance
        %
        if ~isbehav					% is Task PLS
           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
        else						% is Bahavior PLS
            datamatcorrs = rri_corr_maps(behavdata_lst{i}, datamat, n, k);
        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 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];
           datamatcorrs_lst = [datamatcorrs_lst, {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
       [brainlv,s,behavlv] = svd(stacked_datamatcorrs',0);
    else
       [behavlv,s,brainlv] = svd(stacked_datamatcorrs,0);
    end

    original_sal = brainlv * s;
    orig_behavlv = behavlv * s;

    s = diag(s);

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

    % calculate behav scores
    %
    if ~isbehav					% is Task PLS

        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

    else					% is Behavior PLS

        [brainscores, behavscores, lvcorrs] = ...
		rri_get_behavscores(stacked_datamat, stacked_behavdata, ...
		brainlv, behavlv, num_cond_lst(1), num_subj_lst);

        orig_corr = lvcorrs;
        [r1 c1] = size(orig_corr);
        distrib = zeros(r1, c1, num_boot+1);
        distrib(1:r1, 1:c1, 1) = orig_corr;

    end

    rri_progress_ui(progress_hdl,'',1);

    %  Begin bootstrap loop
    %
    %  reorder = rri_mkboot_order(num_cond_lst(1), num_subj_lst, num_boot);

    %  move the following code up, so we don't need to run obserbation
    %  test if bootstrap can not run
    %
    % [reorder, new_num_boot] = ...
	% rri_boot_order(num_subj_lst, num_cond_lst(1), num_boot, 1);

    if new_num_boot ~= num_boot
        num_boot = new_num_boot;
        h0 = findobj(0,'tag','PermutationOptionsFigure');
        h = findobj(h0,'tag','NumBootstrapEdit');
        set(h,'string',num2str(num_boot));

        if isbehav
            distrib = zeros(r1, c1, num_boot+1);
            distrib(1:r1, 1:c1, 1) = orig_corr;
        end
    end

    max_subj_per_group = 8;
    if (sum(num_subj_lst <= max_subj_per_group) == num_groups)
        is_boot_samples = ones(1,num_groups);
    else
        is_boot_samples = zeros(1,num_groups);
    end

    if isempty(reorder)
       return;
    end

    if ~isbehav
        sal_sq = zeros(size(brainlv));
        sal_sum = zeros(size(brainlv));
        dsal_sq = zeros(size(behavlv));
        dsal_sum = zeros(size(behavlv));
    else
        sal_sq = original_sal.^2;
        sal_sum = original_sal;
        dsal_sq = orig_behavlv.^2;
        dsal_sum = orig_behavlv;
    end


    %  Check min% unique values for all behavior variables
    %
    num_LowVariability_behav_boots = zeros(1, size(stacked_behavdata, 2));

    for bw = 1:size(stacked_behavdata, 2)
       for p = 1:num_boot
          v = stacked_behavdata(reorder(:,p),bw);

%          if unique(v) <= size(stacked_behavdata, 1)*0.5
          if rri_islowvariability(v, stacked_behavdata(:,bw))
             num_LowVariability_behav_boots(bw) = num_LowVariability_behav_boots(bw) + 1;
          end
       end
    end

    if any(num_LowVariability_behav_boots)
       disp(' ');
       disp(' ');
       disp('For at least one behavior measure, the minimum unique values of resampled behavior data does not exceed 50% of its total.');
       disp(' ');
       disp(' ');
    end


    for p = 1:num_boot

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

        data_p = stacked_datamat(reorder(:,p),:);
        if 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 ~isbehav
               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);