pls_analysis.m

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

%PLS_ANALYSIS Run PLS analysis on one or more datamats. Datamat list cell
%	array, number of subject list array, and number of condition must
%	be input to the program.
%
%	If this is a Non-Rotated PLS analysis, a contrast text file must
%	be provided;
%
%	If this is a Behavior PLS analysis or Multiblock PLS analysis, 
%	a behavior text file must be provided;
%
%	For behavior text file, columns stand for behavior measures, number
%	of rows for each group is the same as the datamat for that group.
%	All the behavior group should also be stacked together vertically
%	with a total number equal to the sum of the number of rows in
%	datamat list array.
%
%	For contrast text file, columns stand for designs, number of rows
%	for each group is equal to number of conditions. Since all the 
%	groups are stacked vertically, total number of rows in contrast 
%	text file should equal to number of groups multiplied by number 
%	of conditions. For Non-Rotated Behavior PLS, a "condition" in 
%	contrast file corresponds to a condition for a behavior measure.
%	Therefore, for Non-Rotated Behavior PLS, total number of rows in
%	contrast text file should equal to number of groups multiplied by
%	number of conditions multiplied by number of behavior measures.
%
%  Usage: result = pls_analysis(datamat_lst, num_subj_lst, num_cond)
%
%  Here:
%
%	First input parameter 'datamat_lst' is a datamat cell array. Each
%	cell stands for one datamat, which is also referred as a group. 
%	All datamats must be in the form of "subject in condition".
%
%	Input 'num_subj_lst' is an array contains the number of subjects
%	in each group.
%
%	Input 'num_cond' is the number of conditions in datamat_lst.
%
%	Output 'result' is a struct contains all necessary items that are
%	generated based on different scenarios. scenario question will be
%	asked at the beginning of the analysis.
%
%	The 'result' struct could have items like:
%
%	method:			PLS option
%				1. Mean-Centering Task PLS
%				2. Non-Rotated Task PLS
%				3. Regular Behavior PLS
%				4. Multiblock PLS
%				5. Non-Rotated Behavior PLS
%	u:			brainlv or salience
%	s:			singular value
%	v:			designlv or behavlv
%	usc:			brainscores or scalpscores
%	vsc:			designscores or behavscores
%	datamatcorrs_lst:	correlation of behavior data with datamat,
%				only available in behavior PLS.
%	lvcorrs:		correlation of behavior data with usc,
%				only available in behavior PLS.
%	origpost:		posthoc result, only available in behavior
%				PLS, when posthoc file is provided.
%	perm_result:		struct containing permutation result
%		num_perm:	number of permutation
%		sp:		permuted singular value greater than observed
%		sprob:		sp normalized by num_perm
%	boot_result:		struct containing bootstrap result
%		num_boot:	number of bootstrap
%		bootsamp:	bootstrap reorder sample
%		compare_u:	compared salience or compared brain
%		compare_v:	compared designlv or compared behavlv
%		u_se:		standard error of salience or brainlv
%		v_se:		standard error of designlv or behavlv
%		distrib:	orig_usc or orig_corr distribution
%		prop:		orig_usc or orig_corr probability
%
%		following boot_result only available in task PLS:
%
%		orig_usc:	same as usc, with mean-centering on subj
%		ulusc:		upper boundary of orig_usc
%		llusc:		lower boundary of orig_usc
%		ulusc_adj:	percentile of orig_usc distribution with
%				upper boundary of orig_usc
%		llusc_adj:	percentile of orig_usc distribution with
%				lower boundary of orig_usc
%
%		following boot_result only available in behavior PLS:
%
%		orig_corr:	same as lvcorrs
%		ulcorr:		upper boundary of orig_corr
%		llcorr:		lower boundary of orig_corr
%		ulcorr_adj:	percentile of orig_corr distribution with
%				upper boundary of orig_corr
%		llcorr_adj:	percentile of orig_corr distribution with
%				lower boundary of orig_corr
%		num_LowVariability_behav_boots:	display numbers of low
%				variability resampled hehavior data in
%				bootstrap test
%		badbeh:		display bad behav data that is caused by
%				bad re-order (with 0 standard deviation)
%				which will further cause divided by 0
%		countnewtotal:	count the new sample that is re-ordered
%				for badbeh

%   Created on 05-JAN-2005 by Jimmy Shen
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function result = pls_analysis(datamat_lst, num_subj_lst, k)

   result = [];

   %  check input variable
   %
   if nargin ~= 3
      disp(' ');
      disp('Usage: result = pls_analysis(datamat_lst, num_subj_lst, num_cond)');
      disp(' ');disp(' ');
      return;
   end

   if ~iscell(datamat_lst)
      disp(' ');
      disp('Usage: result = pls_analysis(datamat_lst, num_subj_lst, num_cond)');
      disp(' ');
      disp('datamat_lst should be a datamat cell array. Each cell');
      disp('represents a datamat. All datamats must be in the form');
      disp('of "subject in condition".');
      disp(' ');disp(' ');
      return;
   end

   if ~isnumeric(num_subj_lst)
      disp(' ');
      disp('Usage: result = pls_analysis(datamat_lst, num_subj_lst, num_cond)');
      disp(' ');
      disp('num_subj_lst should be a numeric array. Each element');
      disp('represents the number of subjects in the corresponding');
      disp('datamat.');
      disp(' ');disp(' ');
      return;
   end

   if ~isnumeric(k) | length(k) ~= 1
      disp(' ');
      disp('Usage: result = pls_analysis(datamat_lst, num_subj_lst, num_cond)');
      disp(' ');
      disp('num_cond should be a number represents the number of conditions.');
      disp(' ');disp(' ');
      return;
   end

   %  init
   %
   num_groups = length(datamat_lst);
   total_rows = 0;

   %  total rows in stacked datamat (extract datamat from each group,
   %  and stacked them together)
   %
   for g = 1:num_groups
      total_rows = total_rows + size(datamat_lst{g}, 1);
   end


   v7 = version;
   if str2num(v7(1))<7
      singleanalysis = 0;
   else
      singleanalysis = 1;
   end

   pc = computer;
   if singleanalysis & ( strcmp(pc,'GLNXA64') | strcmp(pc,'GLNXI64') | strcmp(pc,'PCWIN64') )

      %  Temporary solution for MATLAB Bug Report ID 268001, 
      %  which has problem to perform single precision math operation 
      %  on any Intel 64-bit machine.
      %
      clc;
      quest = input('\nWe detected that you are running MATLAB on a 64-bit system.\nAccording to MATLAB Bug Report ID 268001, we have to convert\ndata to double precision for Intel based system.\n\nIs this Intel 64-bit machine?\n\n1. No\n2. Yes\n3. Don''t know\n\nPlease select: ','s');
      quest = str2num(quest);
      while isempty(quest) | (quest ~= 1 & quest ~= 2 & quest ~= 3)
         disp(' ');disp(' ');
         disp('Please select either 1, 2, or 3 and then strike ''Enter''');
         quest = input('\nWe detected that you are running MATLAB on a 64-bit system.\nAccording to MATLAB Bug Report ID 268001, we have to convert\ndata to double precision for Intel based system.\n\nIs this Intel 64-bit machine?\n\n1. No\n2. Yes\n3. Don''t know\n\nPlease select: ','s');
         quest = str2num(quest);
      end

      if quest ~= 1
         singleanalysis = 0;
      end

   end;

   for grp=1:length(datamat_lst)
      if singleanalysis
         datamat_lst{grp} = single(datamat_lst{grp});
      else
         datamat_lst{grp} = double(datamat_lst{grp});
      end
   end


   %  which PLS analysis do you choose to run
   %
   clc;
   method = input('\nPLS Option:\n\n1. Mean-Centering Task PLS\n2. Non-Rotated Task PLS\n3. Regular Behavior PLS\n4. Multiblock PLS\n5. Non-Rotated Behavior PLS\n\nPlease select: ','s');
   method = str2num(method);
   while isempty(method) | (method ~= 1 & method ~= 2 & method ~= 3 & method ~= 4 & method ~= 5)
      disp(' ');disp(' ');
      disp('Please select either 1, 2, 3, 4 or 5 and then strike ''Enter''');
      method = input('\nPLS Option:\n\n1. Mean-Centering Task PLS\n2. Non-Rotated Task PLS\n3. Regular Behavior PLS\n4. Multiblock PLS\n5. Non-Rotated Behavior PLS\n\nPlease select: ','s');
      method = str2num(method);
   end

   result.method = method;

   %  which PLS module do you choose to run
   %
   clc;
   disp(' '); disp('Non-Behavior Structure PLS does not permute conditions within a group.');
   module = input('\nPLS Module:\n\n0. Other PLS\n1. Structure PLS\n\nPlease select: ','s');
   module = str2num(module);
   while isempty(module) | (module ~= 0 & module ~= 1)
      disp(' ');disp(' ');
      disp('Please select either 0 or 1 and then strike ''Enter''');
      disp('Non-Behavior Structure PLS does not permute conditions within a group.');
      module = input('\nPLS Module:\n\n0. Other PLS\n1. Structure PLS\n\nPlease select: ','s');
      module = str2num(module);
   end

   is_struct = module;
   result.is_struct = module;

   single_cond_lst = {};

   %  for single condition with multiple groups situation, reconstruct datamat
   %  to make it single group with multiple conditions.
   %
   if method == 1 & k == 1
%      if any(diff(num_subj_lst))
 %        disp('number of subject should be the same for all group in single condition with multiple groups situation.');
  %       disp(' ');disp(' ');
   %      return;
    %  end

%      num_subj_lst = num_subj_lst(1);
 %     k = num_groups;
  %    num_groups = 1;

      tmp = [];

%      for g=1:k
      for g=1:num_groups
         tmp = [tmp; datamat_lst{g}];
      end

      single_cond_lst = {tmp};
   end

   %  For test purpose
   %
   if exist('plslog.m','file')
      switch method
      case 1
         plslog('command-line mean-centering task analysis');
      case 2
         plslog('command-line non-rotated task analysis');
      case 3
         plslog('command-line regular behavior analysis');
      case 4
         plslog('command-line multiblock analysis');
      case 5
         plslog('command-line non-rotated behavior analysis');
      end

      plslog(pwd);
   end

   %  number of permutation
   %
   clc;
   disp(' ');
   num_perm = input('Please enter number of permutation [0]: ','s');
   num_perm = str2num(num_perm);
   if isempty(num_perm) | num_perm <= 0
      num_perm = 0;
   end;

   %  number of bootstrap
   %
   clc;
   disp(' ');
   num_boot = input('Please enter number of bootstrap [0]: ','s');
   num_boot = str2num(num_boot);
   if isempty(num_boot) | num_boot <= 0
      num_boot = 0;
   end;

   %  default value
   %
   contrast_file = '';
   behav_file = '';
   posthoc_file = '';
   stacked_designdata = [];
   stacked_behavdata = [];
   Clim = 95;

   %  need behavior file and posthoc file for regular behav & multiblock PLS
   %
   if method == 3 | method == 4 | method == 5
      clc;
      done = 0;
      while ~done
         disp(' ');
         behav_file = input('Behav filename with path for Behavior PLS: ','s');

         if ~isempty(behav_file) & exist(behav_file,'file')
            try
               stacked_behavdata = load(behav_file);
            catch
               disp(' '); disp('Invalid behavior data file.');
               return;
            end

            if size(stacked_behavdata,1) ~= total_rows
               disp(' '); disp(['Wrong number of rows in behavior data file, which should be ' num2str(total_rows) '.']);
               disp('Please try again.');
            else
               done = 1;
            end
         else
            disp(' '); disp('File does not exist, please try again');
         end
      end

      %  posthoc text file
      %
%      clc;
 %     posthoc_file = input('Posthoc filename with path for Behavior PLS []: ','s');
  %    while ~isempty(posthoc_file) & ~exist(posthoc_file,'file')
   %      disp(' '); disp('Please enter a valid file or strike ''Enter'' if there is no posthoc file');
    %     posthoc_file = input('Posthoc filename with path for Behavior PLS []: ','s');
     % end

      if ~isempty(posthoc_file) & exist(posthoc_file,'file')
         try
            posthoc = load(posthoc_file);

            if size(posthoc,1) ~= size(stacked_behavdata,2)*k*num_groups
               disp(' '); disp('Invalid posthoc row number, posthoc computation will not be performed');
               posthoc_file = '';
               posthoc = [];
               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
   end

   %  confidence level
   %
   clc;
   if num_boot ~= 0
      disp(' ');
      Clim = input('Please enter a confidence level for bootstrap test [95]: ','s');
      Clim = str2num(Clim);
      if isempty(Clim) | Clim < 0 | Clim > 100
         Clim = 95;
      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);