fmri_boot_behav.m

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

function [brainlv,s,behavlv,brainscores,behavscores,lvcorrs, ...
	boot_result,behavdata,evt_list,datamatcorrs_lst]= ...
		fmri_boot_behav(datamat,behavdata,evt_list, ...
		behavdata_lst, newdata_lst, num_subj_lst, ...
		num_boot,num_cond,num_subj,Clim, ...
                min_subj_per_group,is_boot_samples,boot_samples,new_num_boot)

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

    datamatcorrs_lst = {};

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



    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');
%    [reorder, new_num_boot] = rri_boot_order(num_subj_lst, num_cond, num_boot, 0, boot_progress, ...
    [reorder, new_num_boot] = rri_boot_order(num_subj_lst, num_cond, num_boot, 0, ...
	min_subj_per_group,is_boot_samples,boot_samples,new_num_boot);

    if isempty(reorder)
       return;
    end;


    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

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

    s = diag(s);

    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);

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

    rri_progress_ui(progress_hdl,'',1);

    %  Begin bootstrap loop
    %
    % reorder = rri_mkboot_order(num_cond, 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, num_boot);

    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));

        distrib = zeros(r1, c1, num_boot+1);
        distrib(1:r1, 1:c1, 1) = orig_corr;
    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;

    sal_sq = original_sal.^2;
    sal_sum = original_sal;
    dsal_sq = orig_behavlv.^2;
    dsal_sum = orig_behavlv;


    %  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),:);
        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;

            if ~is_boot_samples

		% the code below is mainly trying to find a proper
		% reorder matrix

                % init badbehav cell array to 0
                % which is used to record the bad behav data caused by
                % bad re-order. This var. is for disp only.
                %
                badbehav = zeros(num_cond, size(behavdata_lst{g},2));

                % Check for upcoming NaN and re-sample if necessary.
                % this only happened on behavior analysis, because the
                % 'xcor' inside of '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.
                %
		% For behavpls_boot, also need to account for multiple
		% scans and behavs
		%
                for c=1:k	% traverse all conditions in this group
                    stdmat(c,:) = std(stacked_behavdata(reorder((1+ ...
			(n*(c-1))+span):(n*c+span), p), :));
                end	% scanloop

		% now, check to see if any are zero
		%
                while sum(stdmat(:)==0)>0
                    countnewtotal = countnewtotal + 1;

		    % keep track of scan & behav that force a resample
		    %
                    badbehav(find(stdmat(:)==0)) = ...
			badbehav(find(stdmat(:)==0)) + 1;

                    badbeh{g,countnewtotal} = badbehav;	% save instead of disp

		    % num_boot is just something to be picked to prevent
		    % infinite loop
		    %
                    if countnewtotal > num_boot
%                        msgbox(['countnewtotal exceeds num_boot: ', ...
%                            num2str(countnewtotal)],'modal');
                        disp('Please check behavior data');
                        breakon=1;
                        break;
                    end

                    % reorder(:,p) = rri_mkboot_order(num_cond, ...
                    %    num_subj_lst, 1);
%                    boot_progress = rri_progress_ui('initialize');
%                    reorderp = rri_boot_order(num_subj_lst, num_cond, 1, 0, boot_progress, ...
                    reorderp = rri_boot_order(num_subj_lst, num_cond, 1, 0, ...
			min_subj_per_group,is_boot_samples,boot_samples,1);
                    reorder(:,p) = reorderp;

                    for c=1:k	% recalc stdmat
                        stdmat(c,:) = std(stacked_behavdata(reorder((1+ ...
                            (n*(c-1))+span):(n*c+span), p), :));
                    end	% scanloop

                end	% while

            end		% if ~is_boot_samples

            % now, we can use this proper reorder matrix
            % to generate behav_p & data_p, and then
            % to calculate datamatcoors
            %
            behav_p = stacked_behavdata(reorder(:,p),:);
            data_p = stacked_datamat(reorder(:,p),:);

            if num_groups == 1
               data = rri_corr_maps(behav_p, data_p, n, k);
            else
               data = rri_corr_maps(behav_p(1+span:n*k+span,:), ...
                   data_p(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, sboot, pbehavlv] = svd(stacked_data',0);
        else
           [pbehavlv, sboot, pbrainlv] = svd(stacked_data,0);
        end

        %  rotate pbehavlv to align with the original behavlv
	        %
        rotatemat = rri_bootprocrust(behavlv,pbehavlv);

        %  rescale the bootstrap vectors
        %
        pbrainlv = pbrainlv * sboot * rotatemat;
        pbehavlv = pbehavlv * sboot * rotatemat;

        [brainsctmp, behavsctmp, bcorr] = ...
		rri_get_behavscores(data_p, behav_p, normalize(pbrainlv), ...
		normalize(pbehavlv), num_cond(1), num_subj_lst);

        distrib(1:r1, 1:c1, p+1) = bcorr;

        sal_sum = sal_sum + pbrainlv;
        sal_sq = sal_sq + pbrainlv.^2;

        dsal_sum = dsal_sum + pbehavlv;
        dsal_sq = dsal_sq + pbehavlv.^2;

    end		% for num_boot

    boot_result.num_boot = num_boot;

    boot_result.num_LowVariability_behav_boots = num_LowVariability_behav_boots;

    sal_sum2 = (sal_sum.^2) / (num_boot+1);
    dsal_sum2 = (dsal_sum.^2) / (num_boot+1);

    %  compute standard errors - standard deviation of bootstrap sample
    %  since original sample is part of bootstrap, divide by number of
    %  bootstrap iterations rather than number of bootstraps minus 1
    %
    %  add ceiling to calculations to prevent the following operations
    %  from producing negative/complex numbers
    %
%    brain_se = sqrt((ceil(sal_sq)-ceil(sal_sum2))/(num_boot));
 %   behav_se = sqrt((ceil(dsal_sq)-ceil(dsal_sum2))/(num_boot));

    brain_se = sqrt(((sal_sq)-(sal_sum2))/(num_boot));
    behav_se = sqrt(((dsal_sq)-(dsal_sum2))/(num_boot));

    %  check for zero standard errors - replace with ones
    %
    brain_zeros=find(brain_se<=0);
    if ~isempty(brain_zeros);
        brain_se(brain_zeros)=1;
    end

    behav_zeros=find(behav_se<=0);
    if ~isempty(behav_zeros);
        behav_se(behav_zeros)=1;
    end

    compare = original_sal ./ brain_se;
    compare_behavlv = orig_behavlv ./ behav_se;

    %  for zero standard errors - replace bootstrap ratios with zero
    %  since the ratio makes no sense anyway
    %
    if ~isempty(brain_zeros);
        compare(brain_zeros)=0;
    end

    if ~isempty(behav_zeros);
        compare_behavlv(behav_zeros)=0;
    end

    ul=Clim;
    ll=100-Clim;

        % e.g. 0.05 >> 0.025 for upper & lower tails, two-tailed
        %
        ClimNi = 0.5*(1-(Clim*0.01));

        %  loop to calculate upper and lower CI limits
        %
        for r=1:r1
            for c=1:c1
                ulcorr(r,c)=percentile(distrib(r,c,2:num_boot+1),ul);
                llcorr(r,c)=percentile(distrib(r,c,2:num_boot+1),ll);
                prop(r,c)=length( find(distrib(r,c,2:num_boot+1) ...
                                       <= orig_corr(r,c)) ) / num_boot;
                if prop(r,c)==1 |prop(r,c)==0	% can't calculate the cumulative_gaussian_inv
                    llcorr_adj(r,c)=NaN;
                    ulcorr_adj(r,c)=NaN;
                else

                    % adjusted confidence intervals - in case the
                    % bootstrap samples are extremely skewed

                    % norm inverse to start to adjust conf int
                    %
                    ni=cumulative_gaussian_inv(prop(r,c));

                    % 1st part of recalc the lower conf interval,
                    % this evaluates to +1.96 for 95%CI
                    %
                    uli=(2*ni) + cumulative_gaussian_inv(1-ClimNi);

                    % 1st part of recalc the upper conf interval
                    % e.g -1.96 for 95%CI
                    %
                    lli=(2*ni) + cumulative_gaussian_inv(ClimNi);

                    ncdf_lli=cumulative_gaussian(lli)*100;	% percentile for lower bounds
                    ncdf_uli=cumulative_gaussian(uli)*100;	% percentile for upper bounds

                    % new percentile
                    %
                    llcorr_adj(r,c)=(percentile(distrib(r,c,2:num_boot+1), ...
			ncdf_lli));
                    ulcorr_adj(r,c)=(percentile(distrib(r,c,2:num_boot+1), ...
			ncdf_uli));

                end	% if
            end		% for c
        end		% for r

    boot_result.orig_corr = orig_corr;
    boot_result.ulcorr = ulcorr;
    boot_result.llcorr = llcorr;
    boot_result.ulcorr_adj = ulcorr_adj;
    boot_result.llcorr_adj = llcorr_adj;
    boot_result.prop = prop;
    boot_result.distrib = distrib;
    boot_result.orig_brainlv = original_sal;
    boot_result.brain_se = brain_se;
    boot_result.zero_brain_se = brain_zeros;
    boot_result.compare = compare;
    boot_result.compare_behavlv = compare_behavlv;
    boot_result.bootsamp = reorder;
    boot_result.badbeh = badbeh;
    boot_result.countnewtotal = countnewtotal;

    return;					% fmri_boot_behav