fmri_deviation_perm_test.m

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

function [brainlv,s,designlv,b_scores,d_scores,perm_result,dev_evt_list] = ...
        fmri_deviation_perm_test(st_datamat,num_conditions,evt_list, ...
                                num_perm, subj_group)
%
%  USAGE: [brainlv,s,designlv,perm_result] = ...
%           fmri_deviation_perm_test( st_datamat,num_conditions,evt_list, ...
%				     num_perm, subj_group )
%
%  Apply permutation test on the spatial/temporal datamat
%
%  INPUT:
%    st_datamat -  NxV matrix for the spatial/temporal datamat, where N 
%		   is the number of event onsets and V is the number of 
%		   voxels times the length of temporal window.
%    num_conditions - number of possible conditions
%    evt_list - the event idx for each row in the st_datamat
%    num_perm - number of permutation tests to be performed
%		(set num_perm = 0 for no permutation test)
%    subj_group - a vector, with M elements, contains the number of subjects 
%    		  for the M groups.  
%
%  OUTPUT:
%    brainlv -  a VxC matrix contains the brainlv for each contrast.
%    s -  a CxC diagonal matrix contains the eigenvalue of the 
%         design'*st_datamat matrix.
%    designlv - a CxC matrix contains the designlv for each contrast.    
%    sprod - the ratio of events that exceed the observed s values
%    ** NOTE: For group analysis, the row order for the brainlv and 
%             designlv is corresponding to the condition index in 
%             ascending order from one group to another, ie. 
%             [grp1_cond1 ... grp1_condN grp2_cond1 ...].
%
%   -- Created May 2001 by Wilkin Chau, Rotman Research Institute
%

  progress_hdl = ShowProgress('initialize','Permutation Test Progress');

  msg = sprintf('Working on PLS ... ');
  ShowProgress(progress_hdl,msg)
  if ~isempty(progress_hdl)
     setappdata(progress_hdl,'ProgressScale',1/(num_perm+1)*0.8);
     setappdata(progress_hdl,'ProgressStart',0.2);
     ShowProgress(progress_hdl,0)
  end;


  if isempty(subj_group),
     GROUP_ANALYSIS = 0;			% for nongroup analysis
  else
     GROUP_ANALYSIS = 1;
  end;

  if (GROUP_ANALYSIS == 0),			
%     dev_data = st_datamat - ones(size(st_datamat,1),1)*mean(st_datamat,1);
%     dev_evt_list = evt_list;
     dev_data = gen_dev_data(st_datamat,num_conditions,evt_list);
     dev_evt_list = evt_list;
  else
     dev_data = gen_grp_dev_data(st_datamat,num_conditions,evt_list,subj_group);
%     dev_evt_list = repmat([1:num_conditions],1,length(subj_group));
     dev_evt_list = evt_list;
  end;

  [r c] = size(dev_data);
  if r <= c
     [brainlv,s,designlv] = svd(dev_data',0);
  else
     [designlv,s,brainlv] = svd(dev_data,0);
  end

  s = diag(s);


  if (GROUP_ANALYSIS == 0),			
     d_scores = designlv(evt_list,:);
  else
     d_scores = designlv;

     %  need to expand d_score and reorder with evt_list grp by grp
     %
     d_scores = expand_d_scores(d_scores, num_conditions, evt_list, subj_group);
  end;
  b_scores = st_datamat * brainlv;

  if (num_perm == 0)		% no permutation test
      perm_result = [];
      return;
  end;

  ShowProgress(progress_hdl,1);
  msg = sprintf('Computing permutation orders ...');
  ShowProgress(progress_hdl,msg);

  % generate the permutation orders
  %
  if (GROUP_ANALYSIS == 0),			
     num_repetitions = length(evt_list) / num_conditions;
     perm_order = rri_perm_order(num_repetitions,num_conditions,num_perm);
  else
     perm_order = rri_perm_order(subj_group,num_conditions,num_perm);
  end;

  %---------------------------------------------------------------------
  %-- perform permutation test now -------------------------------------
  %
  sp = zeros(size(s));
  dp = zeros(size(designlv));
  for k=1:num_perm,

    msg = sprintf('Working on permutation ... %d out of %d',k,num_perm);
    ShowProgress(progress_hdl,msg);

    new_order = perm_order(:,k);
    if (GROUP_ANALYSIS == 0),			% for nongroup analysis 
       data_p = gen_dev_data(st_datamat(new_order,:),num_conditions,evt_list);
    else
       data_p = gen_grp_dev_data(st_datamat(new_order,:),num_conditions, ...
                                              evt_list,subj_group);
    end;

    [r c] = size(data_p);
    if r <= c
       [perm_blv,perm_s,perm_dlv] = svd(data_p',0);
    else
       [perm_dlv,perm_s,perm_blv] = svd(data_p,0);
    end

    % rotate designlv to align with the original designlv
    %
    rotatemat = rri_bootprocrust(designlv,perm_dlv);
    perm_dlv = perm_dlv * perm_s * rotatemat;
    perm_s = sqrt(sum(perm_dlv.^2));

    sp = sp+(perm_s'>=s);
    dp = dp+(abs(perm_dlv)>=abs(designlv*(diag(s))));
    
    ShowProgress(progress_hdl,k+1);

  end; 
  %
  %-- permutation loop ------------------------------------------------- 

  sprob = sp./num_perm;
  designlvp = dp./num_perm;

  perm_result.num_perm = num_perm;
  perm_result.permsamp = perm_order;
  perm_result.sp = sp;
  perm_result.s_prob = sprob;
  perm_result.dp = dp;
  perm_result.designlv_prob = designlvp;

  return;


%-------------------------------------------------------------------------
function  data = gen_grp_dev_data(st_datamat,num_conditions,evt_list,subj_group)
%
%  compute the average of st_datamat within the same group
%
  num_group = length(subj_group);
  g_end_idx = 0;
  data = [];

 if num_conditions==1 & num_group>1

   task_idx = cell(1,num_group);

   for g = 1:num_group,
      g_start_idx = g_end_idx + 1;
      g_end_idx = g_start_idx+subj_group(g)-1;
      task_idx{g} = [g_start_idx:g_end_idx];
   end

   %  compute the mean data of each condition for the group 
   %
   if strcmpi(class(st_datamat),'single')
      mean_datamat = single(zeros(num_group,size(st_datamat,2)));
   else
      mean_datamat = zeros(num_group,size(st_datamat,2));
   end

   for i=1:num_group,
      mean_datamat(i,:) =  mean(st_datamat(task_idx{i},:),1);
   end;

   data = mean_datamat - ones(num_group,1)*mean(mean_datamat,1);

 else

  for g = 1:num_group,

     g_start_idx = g_end_idx + 1;
     g_end_idx = g_start_idx+subj_group(g)*num_conditions-1;
     g_range = [g_start_idx:g_end_idx];

     %  store the row indices for each task 
     %
     g_evt_list = zeros(1,length(evt_list));
     g_evt_list(g_range) = evt_list(g_range);

     task_idx = cell(1,num_conditions);
     for i=1:num_conditions,
        task_idx{i} = find(g_evt_list==i);
     end; 

     %  compute the mean data of each condition for the group 
     %
     if strcmpi(class(st_datamat),'single')
        mean_datamat = single(zeros(num_conditions,size(st_datamat,2)));
     else
        mean_datamat = zeros(num_conditions,size(st_datamat,2));
     end

     for i=1:num_conditions,
        mean_datamat(i,:) =  mean(st_datamat(task_idx{i},:),1);
     end;
     grp_datamat = mean_datamat - ones(num_conditions,1)*mean(mean_datamat,1);

     data = [data; grp_datamat];
  end;

 end;

  return; 					% gen_grp_dev_data


%-------------------------------------------------------------------------
function  data = gen_dev_data(st_datamat,num_conditions,evt_list)

  task_idx = cell(1,num_conditions);
  for i=1:num_conditions,
     task_idx{i} = find(evt_list==i);
  end; 

  %  compute the mean data of each condition for the group 
  %
  if strcmpi(class(st_datamat),'single')
      mean_datamat = single(zeros(num_conditions,size(st_datamat,2)));
  else
      mean_datamat = zeros(num_conditions,size(st_datamat,2));
  end

  for i=1:num_conditions,
     mean_datamat(i,:) =  mean(st_datamat(task_idx{i},:),1);
  end;

  data = mean_datamat - ones(num_conditions,1)*mean(st_datamat,1);

  return; 					% gen_dev_data


%-------------------------------------------------------------------------
function hdl = ShowProgress(progress_hdl,info)

  %  'initialize' - return progress handle if any
  %
  if ischar(progress_hdl) & strcmp(lower(progress_hdl),'initialize'),
     if ~isempty(gcf) & isequal(get(gcf,'Tag'),'ProgressFigure'),
         hdl = gcf;
         if ~isempty(info)
             set(hdl,'Name',info);
         end;
     else
         hdl = [];
     end;
     return;
  end;


  if ~isempty(progress_hdl)
     if ischar(info)
         rri_progress_status(progress_hdl,'Show_message',info);
     else
         rri_progress_status(progress_hdl,'Update_bar',info);
     end;
     return;
  end;

  if ischar(info),
%     disp(info)
  end;

  return;					% ShowProgress


%-------------------------------------------------------------------------
function new_d_scores = expand_d_scores(d_scores,num_conditions,evt_list,subj_group)

   new_d_scores = [];
   num_in_grp = [0 num_conditions*subj_group];

   for grp_idx = 1:length(subj_group)
      num_subjs = subj_group(grp_idx);
      first = sum(num_in_grp(1:grp_idx)) + 1;
      last = sum(num_in_grp(1:(grp_idx+1)));
      tmp_evt_list = evt_list(first:last);
      tmp_d_scores = d_scores(((grp_idx-1)*num_conditions+1):grp_idx*num_conditions,:);
      tmp_d_scores = tmp_d_scores(tmp_evt_list,:);	% expand and reorder
      new_d_scores = [new_d_scores;tmp_d_scores];
   end

   return;