fmri_deviation_boot_test.m

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

function [brainlv,s,designlv,b_scores,d_scores,boot_result,dev_evt_list] = ...
        fmri_deviation_boot_test(st_datamat,num_conditions,evt_list, ...
                                num_boot,subj_group, ...
         min_subj_per_group,is_boot_samples,boot_samples,new_num_boot)
%
%  USAGE: [brainlv,s,designlv,boot_result] = ...
%           fmri_deviation_boot_test( st_datamat,num_conditions,evt_list, ...
%				     num_boot, subj_group )
%
%  Apply bootstrap 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_boot - number of bootstraps to be performed
%		(set num_boot = 0 for no bootstrap test)
%    subj_group - a vector, with M elements, contains the number of subjects 
%    		  for one of 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
%

  brainlv = [];
  s = [];
  designlv = [];
  b_scores = [];
  d_scores = [];
  boot_result = [];
  dev_evt_list = [];

  max_subj_per_group = 8;

  if isempty(subj_group),
     num_groups = 0;
     GROUP_ANALYSIS = 0;			% for nongroup analysis

     %  moved from below
     %
     num_rep = length(evt_list) / num_conditions;
 
     % make sure the datamat are blocked by subjects
     [sort_evt,sort_idx] = sort(evt_list);
     new_row_list = reshape(sort_idx,num_rep,num_conditions)';
     new_evt_list = evt_list(new_row_list);

%     boot_progress = rri_progress_ui('initialize');
%     [boot_order, new_num_boot] = rri_boot_order(num_rep,1,num_boot,1,boot_progress, ...
     [boot_order, new_num_boot] = rri_boot_order(num_rep,1,num_boot,1, ...
        min_subj_per_group,is_boot_samples,boot_samples,new_num_boot); 

     if num_rep <= max_subj_per_group
        is_boot_samples = 1;
     else
        is_boot_samples = 0;
     end

  else
     num_groups = length(subj_group);
     GROUP_ANALYSIS = 1;

     %  moved from below
     %
%     boot_progress = rri_progress_ui('initialize');
%     [boot_order, new_num_boot] = rri_boot_order(subj_group,num_conditions,num_boot,1,boot_progress, ...
     [boot_order, new_num_boot] = rri_boot_order(subj_group,num_conditions,num_boot,1, ...
        min_subj_per_group,is_boot_samples,boot_samples,new_num_boot);

     if (sum(subj_group <= length(subj_group)) == num_groups)
        is_boot_samples = 1;
     else
        is_boot_samples = 0;
     end

  end;

  if isempty(boot_order)
     boot_result = [];
     return;
  end

  progress_hdl = ShowProgress('initialize','Bootstrap Progress');

  msg = sprintf('Working on PLS ... ');
  ShowProgress(progress_hdl,msg)
  if ~isempty(progress_hdl)
     setappdata(progress_hdl,'ProgressScale',1/(num_boot+1)*0.8);
     setappdata(progress_hdl,'ProgressStart',0.2);
     ShowProgress(progress_hdl,0)
  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

  original_sal = brainlv*s;
  original_dsal = designlv*s;

  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_boot == 0)		% no bootstrap 
      boot_result = [];
      return;
  end;


  ShowProgress(progress_hdl,1);
  msg = sprintf('Compute bootstrap resampling orders ...');
  ShowProgress(progress_hdl,msg);

%  the following part moved up to prevent observation test to run
%  if bootstrap can not run
%
if(0)

  % generate the resampling orders
  %
  max_subj_per_group = 8;

  if (GROUP_ANALYSIS == 0),			
     num_rep = length(evt_list) / num_conditions;

% following block has been taken care of by rri_boot_order
%     if (num_rep < 3), 
%        boot_result = [];		% must have at least 3 repetitions
%	return; 
%     end;
 
     % make sure the datamat are blocked by subjects
     [sort_evt,sort_idx] = sort(evt_list);
     new_row_list = reshape(sort_idx,num_rep,num_conditions)';
     new_evt_list = evt_list(new_row_list);

%     boot_progress = rri_progress_ui('initialize');
%     [boot_order, new_num_boot] = rri_boot_order(num_rep,1,num_boot,1,boot_progress); 
     [boot_order, new_num_boot] = rri_boot_order(num_rep,1,num_boot,1); 

     if num_rep <= max_subj_per_group
        is_boot_samples = 1;
     else
        is_boot_samples = 0;
     end
  else
%     boot_progress = rri_progress_ui('initialize');
%     [boot_order, new_num_boot] = rri_boot_order(subj_group,num_conditions,num_boot,1,boot_progress);
     [boot_order, new_num_boot] = rri_boot_order(subj_group,num_conditions,num_boot,1);

     if (sum(subj_group <= length(subj_group)) == num_groups)
        is_boot_samples = 1;
     else
        is_boot_samples = 0;
     end
  end;
end

  if new_num_boot ~= num_boot
     num_boot = new_num_boot;
     h0 = findobj(0,'tag','PermutationOptionsFigure');
     h = findobj(h0,'tag','NumBootstrapEdit');
     set(h,'string',num2str(num_boot));

     if ~isempty(progress_hdl)
        setappdata(progress_hdl,'ProgressScale',1/(num_boot+1)*0.8);
     end
  end

  if isempty(boot_order)
     boot_result = [];
     return;
  end

  %---------------------------------------------------------------------
  %-- perform bootstrap now -------------------------------------
  %
  sal_sum = zeros(size(brainlv));
  sal_sq = zeros(size(brainlv));
  dsal_sum = zeros(size(designlv));
  dsal_sq = zeros(size(designlv));

  for k=1:num_boot,

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

    new_order = boot_order(:,k);

    if (GROUP_ANALYSIS == 0),			% for nongroup analysis
       new_row_order = new_row_list(:,new_order);
       data_p = gen_dev_data(st_datamat(new_row_order,:),num_conditions, ...
						  		new_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
       [boot_blv,boot_s,boot_dlv] = svd(data_p',0);
    else
       [boot_dlv,boot_s,boot_blv] = svd(data_p,0);
    end

    % rotate designlv to align with the original designlv
    %
    rotatemat = rri_bootprocrust(designlv,boot_dlv);
    boot_blv = boot_blv * boot_s * rotatemat;
    boot_dlv = boot_dlv * boot_s * rotatemat;

    sal_sum = sal_sum + boot_blv;
    sal_sq = sal_sq +  boot_blv .^ 2;

    dsal_sum = dsal_sum + boot_dlv;
    dsal_sq = dsal_sq + boot_dlv .^ 2;
    
    ShowProgress(progress_hdl,k+1);

  end; 
  %
  %-- bootstrap loop ------------------------------------------------- 


  sal_sum2 = (sal_sum.^2) / num_boot;
  brain_se = sqrt((sal_sq - sal_sum2) / (num_boot - 1));

  dsal_sum2 = (dsal_sum.^2) / num_boot;
  design_se = sqrt((dsal_sq - dsal_sum2) / (num_boot - 1));

  %  check for zero standard errors - replace with ones
  %
  brain_zeros=find(brain_se<=0);
  if ~isempty(brain_zeros);
      brain_se(brain_zeros)=1;
  end

  compare = original_sal ./ brain_se;

  %  for zero standard errors - replace bootstrap ratios with zero
  %  since the ratio makes no sense anyway
  %
  if ~isempty(brain_zeros);
      compare(brain_zeros)=0;
  end

  boot_result.num_boot = num_boot;
  boot_result.bootsamp = boot_order;
  boot_result.brain_se = brain_se;
  boot_result.zero_brain_se = brain_zeros;
  boot_result.design_se = design_se;
  boot_result.compare = compare;
  boot_result.compare_design = original_dsal ./ design_se;

  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, the row order
%  of each group becomes [cond#1 cond#2 ... cond#N]
%
  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;