run_post_process

 run_post_process.m

0001 % run_post_process.m
0002 %
0003 %
0004 
0005 % 2003, Petr Janata
0006 
0007 des_id = 1;
0008 desmat_dir = fullfile(dataroot, '/test/stats/');
0009 
0010 DEMO_STATS = 1;
0011 total_prop_reg  = 0.6;     % The approximate proportion of variance explained by
0012                         % the regressors. It will actually be somewhat more
0013                         % than this
0014 
0015 %get_exp_info;
0016 %opt.title = sprintf('%dSil, TR=%1.1f, note_dur=%d ms, cue_stim_soa=%d-%d s, ITI=%1.1f-%1.1f s, total_dur=%2.1f min',tinfo.num_trials(ntt),TR/1000,ms_per_note, min(cue_stim_soa_range)/1000, max(cue_stim_soa_range)/1000, min(inter_trial_interval_range)/1000, max(inter_trial_interval_range)/1000, total_time_min);
0017 
0018 %opt.title = strrep(opt.title,'_','\_');
0019 opt.printfig = 0;
0020 opt.append_str = '-append';
0021 %opt.append_str = '';
0022 opt.figstub = fullfile(desmat_dir,sprintf('pb_des%02d', des_id));
0023 
0024 % Create a plot that shows the correlations among the regressors
0025 [corrmat] = post_process_corrmat(fullfile(desmat_dir,'SPM_fMRIDesMtx'),opt);
0026 
0027 if DEMO_STATS
0028   
0029   % Load the design matrix that we created
0030   load(fullfile(desmat_dir,'SPM_fMRIDesMtx'))
0031   X = xX.X;
0032   
0033   % Remove the constant term that was added by SPM99
0034   X(:,end) = [];
0035   
0036   % Create some data.
0037   % Assume that each of the regressors explains an equal amount of variance.
0038   
0039   nreg = size(X,2);
0040   prop_per_reg = total_prop_reg/nreg
0041   
0042   % In creating the fake data, weight some of the regressors of we want to
0043   weights = ones(1,nreg);
0044   weights(2:2:end) = 1; % use -1 to make some of the regressors negatively correlated
0045   
0046   fakeX = X .* repmat(weights,size(X,1),1);
0047   
0048   % Add the regressors together and scale them to between 0 and 1.  This is the
0049   % part of the fake data than can be explained by our model
0050   reg_sum = sum(fakeX,2);
0051   reg_sum = reg_sum-min(reg_sum);
0052   reg_sum = reg_sum/max(reg_sum);
0053   
0054   % Create a noise vector of values between 0 and 1.  This is the part of the
0055   % data that can't be explained by the model (beyond some chance level)
0056   rand_vect = rand(size(reg_sum));
0057   
0058   % Add the signal and noise together in desired proportion
0059   fake_data = reg_sum*sqrt(total_prop_reg) +rand_vect*sqrt(1-total_prop_reg);
0060 
0061   % Create a plot that shows the sum of the regressors and the fake data
0062   figure(3), clf
0063   subplot(2,1,1)
0064   plot(reg_sum)
0065   
0066   subplot(2,1,2)
0067   plot(fake_data) 
0068 end
0069 
0070 origX = X;
0071 
0072 %
0073 % Do some generalized linear model fitting
0074 %
0075 
0076 % Fit only a subset of the regressors
0077 use_reg = [1];
0078 X = origX(:,use_reg);
0079 [b_sub, dev, stats_sub] = glmfit(zscore(X), zscore(fake_data), 'normal');
0080 fprintf('Using regressors: %s\n', sprintf('%d ',use_reg));
0081 fprintf('Proportion of variance explained: %s\n', sprintf('%1.2e, ', b_sub(2:end).^2));
0082 fprintf('Significance (p<): %s\n', sprintf('%1.4f, ', stats_sub.p(2:end)))
0083 
0084 % Calculate the proportion of variance explained
0085 
0086 % Fit the whole model
0087 X = origX;
0088 [b, dev, stats] = glmfit(zscore(X), zscore(fake_data), 'normal');
0089 
0090 z = zscore(fake_data);
0091 tss = z'*z;
0092 ess = stats.resid'*stats.resid;
0093 
0094 fprintf('\nUsing all regressors\n');
0095 fprintf('Total SS via beta^2: %1.3f\n', sum(b(2:end).^2));
0096 fprintf('Total SS via resid: %1.3f\n', 1-ess/tss);
0097 fprintf('Proportion of variance explained: %s\n', sprintf('%1.2e, ', b(2:end).^2));
0098 fprintf('Significance (p<): %s\n', sprintf('%1.4f, ', stats.p(2:end)))
0099 
0100 % Create a plot
0101 bmat = zeros(length(b)-1,2);
0102 pmat = zeros(size(bmat));
0103 
0104 bmat(use_reg,1) = b_sub(2:end);
0105 bmat(:,2) = b(2:end);
0106 
0107 pmat(use_reg,1) = stats_sub.p(2:end);
0108 pmat(:,2) = stats.p(2:end);
0109 
0110 figure(4),clf
0111 subplot(2,1,1)
0112 bar(bmat)
0113 title('Beta values')
0114 
0115 subplot(2,1,2)
0116 bar(pmat)
0117 l = line([get(gca,'xlim')],[0.05 0.05]); % draw a red line at a p<0.05 threshold
0118 set(l,'color','r')
0119 title('P values')