pls_analysis.m

matlab里的PLS回归分析,含英文说明,比较复杂

               disp(' '); disp('Press any key to continue');
               pause;
            end
         catch
            disp(' '); disp('Invalid posthoc data file, posthoc computation will not be performed');
            posthoc_file = '';
            posthoc = [];
            disp(' '); disp('Press any key to continue');
            pause;
         end
      end

      %  confidence level
      %
      clc;
      if num_boot ~= 0
         Clim = input('Please enter a confidence level for Behavior PLS [95]: ','s');
         Clim = str2num(Clim);
         if isempty(Clim) | Clim < 0 | Clim > 100
            Clim = 95;
         end;
      end
   end

   %  need condition deselect for multiblock PLS
   %
   if method == 4
      k_str = num2str(1:k);
      msg='Please select index of only those conditions for behavior block: ';
      tmp=inputdlg({msg}, 'Behavior block conditions deselect', 1, {k_str});

      if isempty(tmp)
         bscan = 1:k;
      else
         bscan = round(str2num(tmp{1}));

         if isempty(bscan) | max(bscan) < 1 | max(bscan) > k
            bscan = 1:k;
         end
      end

      result.bscan = bscan;
      kk = length(bscan);

      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

   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);
      else
         [min_subj_per_group, is_boot_samples, boot_samples, new_num_boot] ...
		= rri_boot_check(num_subj_lst, k, num_boot);
      end
   end

   %  init variable for the following loop
   %
   vsc = [];
   datamatcorrs_lst = {};
   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

         TBdatamatcorrs = [];
      case 2
         datamatcorrs = rri_task_mean(datamat,n);
         TBdatamatcorrs = [];
      case 3
         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];
      end

   end	% for num_groups

   fprintf('\n');

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

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

   if method == 2	% different computation for non-rotate task

      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, 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}
      usc = stacked_datamat * u;
      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

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