1 function example_djm_unbalanced()
   2 %{
   3 %%%%%%%%%%%%%%%%%%%%%
   4 djm: This shows a simple group analysis for the vt mask of the AK6 dataset.
   5 It illustrates the issue with unbalanced samples, and three possible solutions.
   6 (Note this code could be made more efficient be rearranging the loops, but
   7 instead I've tried to keep it reasonably simple.)
   8 %%%%%%%%%%%%%%%%%%%%%
   9 %}
  10 
  11 %% Define data
  12 config=cosmo_config();
  13 subs={'s01','s02','s03','s04','s05','s06','s07','s08'};
  14 outputs={'accuracy','balanced_accuracy','acc_rebal'};
  15 % The first two of these are possible arguments to cosmo_crossvalidation_measure
  16 % The last one will use "accuracy" after calling "cosmo_balance_partitions"
  17 
  18 figure(10); clf(10);
  19 
  20 rng(0); % reset random number generator for reproducability
  21 
  22 for balanced=[true false]
  23   for signal=[true false]
  24     situation=sprintf('Signal=%g; balanced=%g',signal,balanced);
  25     fprintf('\n%s',situation) % display current situation to command line
  26 
  27     result=nan( numel(subs), numel(outputs) ); % djm, preallocate array
  28     for s=1:numel(subs)
  29       fprintf('.')
  30       data_path=fullfile(config.tutorial_data_path,'ak6',subs{s});
  31       data_fn=fullfile(data_path,'glm_T_stats_perrun.nii');
  32       mask_fn=fullfile(data_path,'vt_mask.nii');
  33       ds=cosmo_fmri_dataset(data_fn,'mask',mask_fn,...
  34       'targets',repmat(1:6,1,10),... % 6 types of animal
  35       'chunks',floor(((1:60)-1)/6)+1); % 10 runs
  36 
  37       % remove constant features (due to liberal masking)
  38       ds=cosmo_remove_useless_data(ds);
  39 
  40       % Assign a function handle to the cosmo_crossvalidation_measure
  41       measure=@cosmo_crossvalidation_measure;
  42 
  43       % Make a struct containing the arguments for the measure:
  44       args=struct();
  45       args.classifier=@cosmo_classify_naive_bayes; % use naive_bayes in this example because lda insists on balanced samples
  46       args.partitions=cosmo_nfold_partitioner(ds);
  47 
  48       % if signal==false, replace real data with random data
  49       if ~signal, ds.samples=randn(size(ds.samples)); end
  50 
  51       % reassign targets into two classes, and assume chance is 0.5:
  52       if balanced,
  53         ds.sa.targets=  (ds.sa.targets<4)+1; % larger versus smaller animals
  54       else
  55         ds.sa.targets=  (ds.sa.targets<3)+1; % primates vs others
  56         args.check_partitions=false; % otherwise CoSMoMVPA should spot this and return a helpful error message
  57       endif
  58 
  59       for o=1:numel(outputs)
  60 
  61         if strcmp(outputs{o},'acc_rebal')
  62           args.output='accuracy';
  63           args.partitions=cosmo_balance_partitions(args.partitions,ds);
  64         else
  65           args.output=outputs{o};
  66         endif
  67 
  68         % Apply the measure to ds, with args as second argument.
  69         ds_accuracy=measure(ds,args);
  70 
  71         % store the result as accuracy minus chance (so positve values indicate information present)
  72         result(s,o)=ds_accuracy.samples-0.5;
  73 
  74       end % next output score
  75     end % next subject
  76 
  77     % Run ttest versus chance (now zero)
  78     [h, p, ci, stat]=ttest(result);
  79 
  80     % plot results with 95% confidence intervals
  81     figure(10); subplot(2,2,signal*2+balanced+1)
  82     errorbar(mean(result),range(ci)./2,'o')
  83     set(gca,'xtick',1:numel(outputs),'xticklabel',outputs)
  84     ylabel('accuracy minus chance')
  85     title(situation)
  86     ylim([-0.2 0.5]);
  87     hold on
  88     plot(xlim,[0 0],'--')
  89     drawnow
  90 
  91     % pause to discuss
  92     %keyboard
  93 
  94   end % signal present?
  95 end %balanced samples?
  96 
  97 return