pls_analysis.m

国外的一个PLStoolbox,主要用于处理图象,也可以用来回归,欢迎使用

            data_p = stacked_datamat(Treorder(:,p),:);
            behav_p = stacked_behavdata(Breorder(:,p),:);
         elseif method == 3 | method == 5
            behav_p = stacked_behavdata(reorder(:,p),:);
         else
            data_p = stacked_datamat(reorder(:,p),:);
         end

         stacked_TBdata_p = [];
         stacked_data_p = [];

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

            switch method
            case 1
               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);	% pet, erp & multiblock style
                  else
%                     meanmat = rri_task_mean(data_p(1+span:n*k+span,:),n);
%                     data = meanmat-(ones(k,1)*mean(meanmat));
                     data = 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
               elseif g==1
                  data = rri_task_mean1(data_p,num_subj_lst)-ones(num_groups,1)*mean(data_p);
               end

               TBdata = [];
            case 2
               if num_groups == 1
                  data = rri_task_mean(data_p, n);
               else
                  data = rri_task_mean(data_p(1+span:n*k+span,:), n);
               end

               TBdata = [];
            case {3, 5}
		% 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'];

                       disp(' '); disp(msg);
                       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

                TBdata = [];
            case 4
               if num_groups == 1
%                  meanmat = rri_task_mean(data_p,n);
%                  Tdata = meanmat-(ones(k,1)*mean(meanmat));
                  Tdata = rri_task_mean(data_p,n)-ones(k,1)*mean(data_p);	% pet, erp & multiblock style
               else
%                  meanmat = rri_task_mean(data_p(1+span:n*k+span,:),n);
%                  Tdata = meanmat-(ones(k,1)*mean(meanmat));
                  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

                %  no comments below. because it is the same as case 3 above
                %
                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'];

                       disp(' '); disp(msg);
                       return;
                    end
                end

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

            if isempty(single_cond_lst) | g==1
               stacked_TBdata_p = [stacked_TBdata_p; TBdata];
               stacked_data_p = [stacked_data_p; data];	% normalized TB
            end
         end

         if method == 2 | method == 5
            crossblock = normalize(stacked_designdata)'*stacked_data_p;
            sperm = sqrt(sum(crossblock.^2,2));
            sp = sp + (sperm >= s);
         else
            %  Singular Value Decomposition
            %
            [r c] = size(stacked_data_p);
            if r <= c
               [pu, sperm, pv] = svd(stacked_data_p',0);
            else
               [pv, sperm, pu] = svd(stacked_data_p,0);
            end

            %  rotate pv to align with the original v
            %
            rotatemat = rri_bootprocrust(v, pv);

            %  rescale the vectors
            %
            pv = pv * sperm * rotatemat;

            sperm = sqrt(sum(pv.^2));

            if method ~= 4
               sp = sp + (sperm'>=s);

               if ~isempty(posthoc_file)
                  tmp = rri_xcor(posthoc, pv);
                  porigpost = porigpost + (abs(tmp) >= abs(origpost));
               end
            else
               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_file)
                  tmp = rri_xcor(posthoc, Bpv);
                  porigpost = porigpost + (abs(tmp) >= abs(origpost));
               end
            end
         end		% if method

         if method == 4
            ptotal_s = sum(stacked_TBdata_p(:).^2);
            per = diag(sperm).^2 / sum(diag(sperm).^2);
            sperm = sqrt(per * ptotal_s);
            pv = normalize(pv) * diag(sperm);
            sp = sp + (sperm>=org_s);
            vp = vp + (abs(pv) >= abs(org_v));
         end		% if method 4

         if pcntacc > 70
            fprintf('\n');
            pcntacc = 0;
         end

      end		% for num_perm

      fprintf('\n');

      if method == 4
         result.perm_result.permsampT = Treorder;
         result.perm_result.permsampB = Breorder;
      else
         result.perm_result.permsamp = reorder;
      end

      result.perm_result.num_perm = num_perm;
      result.perm_result.sp = sp;
      result.perm_result.sprob = sp ./ (num_perm + 1);

      if method == 4
         result.perm_result.vprob = vp ./ (num_perm + 1);
      end

   end	% if perm

   %-------------------------_______________________-----------------------

   if num_boot > 0

      %  keeps track of number of times a new bootstrap had to be generated
      %
      countnewtotal=0;
      num_LowVariability_behav_boots = [];
      badbeh = [];
      fprintf('\nMaking resampling matrix for bootstrap ...\n');

      %  include original un-resampled order only for mean-centering task PLS
      %
      if method == 1
         [reorder, new_num_boot] = rri_boot_order(num_subj_lst,k,num_boot, ...
            1,min_subj_per_group,is_boot_samples,boot_samples,new_num_boot,1);
      else
         [reorder, new_num_boot]=rri_boot_order(num_subj_lst,k,num_boot, ...
            0,min_subj_per_group,is_boot_samples,boot_samples,new_num_boot,1);
      end

      if isempty(reorder)
         disp('bootstrap order is not available');
         return;
      end;

      if new_num_boot ~= num_boot
         num_boot = new_num_boot;
      end

      if method == 3 | method == 4 | method == 5
         orig_corr = lvcorrs;

         [r1 c1] = size(orig_corr);
         distrib = zeros(r1, c1, num_boot+1);
         distrib(:, :, 1) = orig_corr;
      elseif method == 1 | method == 2
         orig_usc = [];

         first = 1;
         last = 0;

         for i = 1:length(num_subj_lst)
            last = last + k*num_subj_lst(i);
            orig_usc = [orig_usc; rri_task_mean(usc2(first:last,:), num_subj_lst(i))];
            first = last + 1;
         end

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

      max_subj_per_group = 8;
      if (sum(num_subj_lst <= max_subj_per_group) == num_groups)
         is_boot_samples = 1;
      else
         is_boot_samples = 0;
      end

      switch method
      case 1
         u_sq = zeros(size(u));
         u_sum = zeros(size(u));
         v_sq = zeros(size(v));
         v_sum = zeros(size(v));
      case {2, 5}
         u_sq = u.^2;
         u_sum = u;
         v_sq = v.^2;
         v_sum = v;
      case {3, 4}
         u_sq = original_u.^2;
         u_sum = original_u;
         v_sq = original_v.^2;
         v_sum = original_v;
      end

      if method==3 | method==4 | method==5
         %  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
               vv = stacked_behavdata(reorder(:,p),bw);

%               if unique(vv) <= size(stacked_behavdata, 1)*0.5
               if rri_islowvariability(vv, 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
      end

      pcntacc = fprintf('Working on %d bootstraps:', num_boot);

      for p=1:num_boot

         pcntacc = pcntacc + fprintf(' %d', p);

         data_p = stacked_datamat(reorder(:,p),:);
         stacked_data_p = [];
         stacked_smeanmat_p = [];

         for g=1:num_groups
            n = num_subj_lst(g);
            span = sum(num_subj_lst(1:g-1)) * k;	% group length

            switch method
            case 1
               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);	% pet, erp & multiblock style
                     smeanmat_p = data_p-ones(size(data_p,1),1)*mean(data_p);
                  else
%                     meanmat = rri_task_mean(data_p(1+span:n*k+span,:),n);
%                     data = meanmat-(ones(k,1)*mean(meanmat));
                     data = 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
                     smeanmat_p = data_p(1+span:n*k+span,:)-ones(size(data_p(1+span:n*k+span,:),1),1)*mean(data_p(1+span:n*k+span,:));
                  end
               elseif g==1
                  data = rri_task_mean1(data_p,num_subj_lst)-ones(num_groups,1)*mean(data_p);
                  smeanmat_p = data_p-ones(size(data_p,1),1)*mean(data_p);
               end
            case 2
               if num_groups == 1
                  smeanmat_p = data_p-ones(size(data_p,1),1)*mean(data_p);	% pet, erp & multiblock style
                  data = rri_task_mean(data_p, n);
               else
                  smeanmat_p = data_p(1+span:n*k+span,:)-ones(size(data_p(1+span:n*k+span,:),1),1)*mean(data_p(1+span:n*k+span,:));	% pet, erp & multiblock style
                  data = rri_task_mean(data_p(1+span:n*k+span,:), n);
               end
            case {3, 5}
                if ~is_boot_samples

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

                   % init badbehav array to 0
                   % which is used to record the bad behav data caused by
                   % bad re-order. This variable is for disp only.
                   %
                   badbehav = zeros(k, size(stacked_behavdata,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
                           disp('Please check behavior data');
                           breakon=1;
                           break;
                       end

                       reorderp = rri_boot_order(num_subj_lst,k,1, ...
                           0,min_subj_per_group,is_boot_samples, ...