pls_analysis.m

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

      if ~isempty(posthoc_file)
         row_idx = [];

         for g = 1:num_groups
            tmp = reshape(1:size(stacked_behavdata,2)*k, [size(stacked_behavdata,2) k]);
            tmp = tmp(:,bscan);
            row_idx = [row_idx tmp(:)+size(stacked_behavdata,2)*k*(g-1)];
         end

         row_idx = row_idx(:);
         posthoc = posthoc(row_idx,:);
      end
   end

   %  need contrast file for non-rotated task & non-rotated behav PLS
   %
   if method == 2 | method == 5
      clc;
      done = 0;
      while ~done
         disp(' ');
         contrast_file = input('Contrast filename with path for Non-rotated PLS: ','s');

         if ~isempty(contrast_file) & exist(contrast_file,'file')
            try
               stacked_designdata = load(contrast_file);
            catch
               disp(' '); disp('Invalid contrast data file.');
               return;
            end

            if method == 2 & size(stacked_designdata,1) ~= num_groups * k
               disp(' '); disp(['Wrong number of rows in contrast data file, which should be ' num2str(num_groups * k) '.']);
               disp('Please try again.');
            elseif method == 5 & size(stacked_designdata,1) ~= num_groups * k * size(stacked_behavdata,2)
               disp(' '); disp(['Wrong number of rows in contrast data file, which should be ' num2str(num_groups * k * size(stacked_behavdata,2)) '.']);
               disp('Please try again.');
            else
               done = 1;
            end
         else
            disp(' '); disp('File does not exist, please try again');
         end
      end

      %  check if the contrast matrix is rank deficient
      %
      if (rank(stacked_designdata) ~= size(stacked_designdata,2))
         disp('Your Contrast matrix is rank deficient.');
      end;

      %  check if the contrast matrix is orthogonal
      %
      check_orth = abs(triu(stacked_designdata'*stacked_designdata) - tril(stacked_designdata'*stacked_designdata));
      if max(check_orth(:)) > 1e-4
         disp('Effects expressed by each contrast are not independent. Check variable');
         disp('lvintercorrs in result file to see overlap of effects between LVs');
      end
   end

   if num_boot > 0
      if method == 1
         [min_subj_per_group, is_boot_samples, boot_samples, new_num_boot] ...
		= rri_boot_check(num_subj_lst, k, num_boot, 1, 0, 1);
      else
         [min_subj_per_group, is_boot_samples, boot_samples, new_num_boot] ...
		= rri_boot_check(num_subj_lst, k, num_boot, 0, 0, 1);
      end
   end

   %  init variable for the following loop
   %
   vsc = [];
   datamatcorrs_lst = {};
   stacked_smeanmat = [];
   stacked_datamat = [];
   stacked_datamatcorrs = [];
   stacked_TBdatamatcorrs = [];	% multiblock PLS before normalize
   first = 1;			% first row for stacked_behavdata

   %  loop accross the groups, and
   %  calculate datamatcorrs for each group
   %
   fprintf('Stacking datamat from group:');

   for g = 1:num_groups

      fprintf(' %d', g);

      if isempty(single_cond_lst)
         datamat = datamat_lst{g};
         stacked_datamat = [stacked_datamat; datamat];
      elseif g==1
         datamat = single_cond_lst{1};
         stacked_datamat = [stacked_datamat; datamat];	% only do once
      end

      n = num_subj_lst(g);

      % compute correlation or covariance
      %
      switch method
      case 1
%         meanmat = rri_task_mean(datamat,n);
%         datamatcorrs = meanmat-(ones(k,1)*mean(meanmat));

         if isempty(single_cond_lst)
            datamatcorrs = rri_task_mean(datamat,n)-ones(k,1)*mean(datamat);	% pet, erp & multiblock style
         elseif g==1
            datamatcorrs = rri_task_mean1(datamat,num_subj_lst)-ones(num_groups,1)*mean(datamat);
         end

         smeanmat = datamat-ones(size(datamat,1),1)*mean(datamat);
         TBdatamatcorrs = [];
      case 2
         if num_boot > 0
            smeanmat = datamat-ones(size(datamat,1),1)*mean(datamat);	% calc orig_usc for boot CI
         end

         datamatcorrs = rri_task_mean(datamat,n);
         TBdatamatcorrs = [];
      case {3, 5}
         last = first + size(datamat,1) - 1;
         behavdata = stacked_behavdata(first:last, :);
         first = last + 1;
         datamatcorrs = rri_corr_maps(behavdata, datamat, n, k);
         datamatcorrs_lst = [datamatcorrs_lst, {datamatcorrs}];
         TBdatamatcorrs = [];
      case 4
%         meanmat = rri_task_mean(datamat,n);
%         Tdatamatcorrs = meanmat-(ones(k,1)*mean(meanmat));
         Tdatamatcorrs = rri_task_mean(datamat,n)-ones(k,1)*mean(datamat);	% pet, erp & multiblock style

         last = first + size(datamat,1) - 1;
         behavdata = stacked_behavdata(first:last, :);
         first = last + 1;
         Bdatamatcorrs = rri_corr_maps_notall(behavdata, datamat, n, bscan);
         datamatcorrs_lst = [datamatcorrs_lst, {Bdatamatcorrs}];

         %  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 differences
         %
         datamatcorrs = [normalize(Tdatamatcorrs,2); normalize(Bdatamatcorrs,2)];
      end

      if isempty(single_cond_lst) | g==1
         stacked_TBdatamatcorrs = [stacked_TBdatamatcorrs; TBdatamatcorrs];
         stacked_datamatcorrs = [stacked_datamatcorrs; datamatcorrs];

         if num_boot > 0 & (method == 1 | method == 2)
            stacked_smeanmat = [stacked_smeanmat; smeanmat];
         end
      end

   end	% for num_groups

   fprintf('\n');

   %  save datamatcorrs_lst, if behav PLS
   %
   if method == 3 | method == 4 | method == 5
      result.datamatcorrs_lst = datamatcorrs_lst;
   end

   clc;
   disp(' '); disp('Calculating LVs ...');

   if method==2 | method==5	% different computation for non-rotated pls

      crossblock = normalize(stacked_designdata)'*stacked_datamatcorrs;
      u = crossblock';
      s = sqrt(sum(crossblock.^2, 2));
      v = stacked_designdata;
      normalized_u = normalize(u);
      lvintercorrs = normalized_u'*normalized_u;
      result.lvintercorrs = lvintercorrs;

   else				% mean-centering, regular behavior or multiblock

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

      s = diag(s);

   end

   org_s = s;
   org_v = v;

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

   %  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)
   %
   if method == 4

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

   end

   %  save u,s,v for all situations (observation)
   %
   result.u = u;
   result.s = s;
   result.v = v;

%   clc;
   disp('Calculating Scores ...');

   switch method
   case {1, 2}
      if method == 1
         usc = stacked_datamat * u;

         if num_boot > 0
            usc2 = stacked_smeanmat * u;
         end
      elseif method == 2
         usc = stacked_datamat * normalized_u;

         if num_boot > 0
            usc2 = stacked_smeanmat * normalized_u;
         end
      end

      num_col = size(v, 2);

      % expand the num_subj for each row (cond)
      % did the samething as testvec
      %
      for g = 1:num_groups
         n = num_subj_lst(g);

         tmp = reshape(v((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]);

         vsc = [vsc; tmp];			% stack by groups
      end
   case 3
      [usc, vsc, lvcorrs] = rri_get_behavscores(stacked_datamat, ...
		stacked_behavdata, u, v, k, num_subj_lst);

      result.lvcorrs = lvcorrs;

      if ~isempty(posthoc_file)
         origpost = rri_xcor(posthoc, v);
         porigpost = zeros(size(origpost));
         result.origpost = origpost;
      end
   case 4

      %  Separate v into 2 parts: Tv (Task) and Bv (Behavior)
      %
      Tv = [];
      Bv = [];

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

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

      num_col = size(Tv, 2);

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

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

         tmp = reshape(Tv((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]);

         Tvsc = [Tvsc; 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);
         row_idx = [row_idx ; tmp(:) + last];
         last = last + n*k;
      end

      Tusc = stacked_datamat * u;

      [Busc, Bvsc, lvcorrs] = rri_get_behavscores(stacked_datamat(row_idx,:), ...
		stacked_behavdata(row_idx,:), u, Bv, kk, num_subj_lst);

      usc = [Tusc; Busc];
      vsc = [Tvsc; Bvsc];
      result.TBusc = {Tusc, Busc};
      result.TBvsc = {Tvsc, Bvsc};
      result.TBv = {Tv, Bv};

%      lvcorrs = original_v;    (?)
      result.lvcorrs = lvcorrs;

      if ~isempty(posthoc_file)
         origpost = rri_xcor(posthoc, Bv);
         porigpost = zeros(size(origpost));
         result.origpost = origpost;
      end

   case 5
      [usc, vsc, lvcorrs] = rri_get_behavscores(stacked_datamat, ...
		stacked_behavdata, normalized_u, v, k, num_subj_lst);

      result.lvcorrs = lvcorrs;

      if ~isempty(posthoc_file)
         origpost = rri_xcor(posthoc, v);
         porigpost = zeros(size(origpost));
         result.origpost = origpost;
      end
   end

   %  save Scores for all situations
   %
   result.usc = usc;
   result.vsc = vsc;

   %  save Input for all situations
   %
   if ~isempty(stacked_behavdata)
      result.behavdata = stacked_behavdata;
   end

   if ~isempty(stacked_designdata)
      result.designdata = stacked_designdata;
   end

   result.datamat_lst = datamat_lst;
   result.num_subj_lst = num_subj_lst;
   result.num_conditions = k;

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

   if num_perm > 0

      sp = zeros(size(s));
      vp = zeros(size(v));
      fprintf('\nMaking resampling matrix for permutation ...\n');

      if method == 3 | method == 5
         for p = 1:num_perm
            reorder(:,p) = [randperm(size(stacked_datamat,1))'];
         end
      elseif method == 4
         Treorder = rri_perm_order(num_subj_lst, k, num_perm, is_struct);

         for p = 1:num_perm
            Breorder(:,p) = [randperm(size(stacked_datamat,1))'];
         end
      else
         reorder = rri_perm_order(num_subj_lst, k, num_perm, is_struct);
      end

      pcntacc = fprintf('Working on %d permutations:', num_perm);

      for p = 1:num_perm

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

         if method == 4