pet_multiblock_boot.m

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

        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_lst(g);
            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);	% pet, erp & multiblock style
            else
               Tdata = rri_task_mean(data_p(1+span:n*k+span,:),n)-ones(k,1)*mean(data_p(1+span:n*k+span,:));	% pet, erp & multiblock style
            end

            if 1

                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_lst(1), ...
                       %    num_subj_lst, 1);
%                       boot_progress = rri_progress_ui('initialize');
%                       reorderp = rri_boot_order(num_subj_lst, num_cond_lst(1), 1, 0, boot_progress, ...
                       reorderp = rri_boot_order(num_subj_lst, num_cond_lst(1), 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
                   Bdata = rri_corr_maps_notall(behav_p, data_p, n, bscan);
                else
                   Bdata = rri_corr_maps_notall(behav_p(1+span:n*k+span,:), ...     
				data_p(1+span:n*k+span,:), n, bscan);
                end
            end

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

            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(v,pbehavlv);

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

            [brainsctmp, behavsctmp, bcorr] = ...
		rri_get_behavscores(data_p(row_idx, :), ...
		behav_p(row_idx, :), ...
		normalize(pbrainlv), ...
		normalize(pbehavlv), kk, num_subj_lst);

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

        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.bootsamp = reorder;

    boot_result.num_LowVariability_behav_boots = num_LowVariability_behav_boots;

    if 1

        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(((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.badbeh = badbeh;
        boot_result.countnewtotal = countnewtotal;

    end		% enf if 1

    return;					% pet_multiblock_boot


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   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