Home > functions > job > vb_job_chi2test_pop1.m

vb_job_chi2test_pop1

PURPOSE ^

Chi2 test for testing significance of two time-windows of one datasets

SYNOPSIS ^

function vb_job_chi2test_pop1(test_parm);

DESCRIPTION ^

 Chi2 test for testing significance of two time-windows of one datasets

 function vb_job_chi2test_pop1(test_parm);

 Field of 'test_parm'
 ---- Input
 currfile : result of VB estimation. 
 tsmplpre  : pre-trigger time points (time origin is 0)
             (e.g. [] (1 sample chi2 test), [1:100]')
 tsmplpost : post-trigger time points (time origin is 0)
             (e.g. [1:100]', [1:300;]')

 ---- Optional Input  (default value)
 testsize (0.05) : Significance level before Bonferroni correction 
                   (e.g. 0.05, 0.01 and so on)
 multiple_correction (ON) : Bonferroni correction is applied if ON
 analyzefile ([]): MRI analyze file. If exists, the result will be
                   displayed on MRI transvese images
 version   (['ANALYZE'])   : Select 'ANALYZE' or 'SBI' 
                   if your vbemge version is ver0.2* or later, choose 'ANALYZE',       
                   otherwise 'SBI'
 outputfile ([])  : If specified, the variables 'pval','ix_dip','test_parm' are saved.  
                      

 --- Example
 test_parm.currfile = ['c00015a-a_UR4.curr.mat'];
 test_parm.tsmplpre = [1:300]'; %pre-trigger
 test_parm.tsmplpost = [360:420]'; %post-trigger 
 test_parm.analyzefile = '/home/dcfs2/oyamashi/matlab/vbmeg2/realdata/retino/SH/data/3D.img';
 test_parm.version = 'ANALYZE'

 ********** Notice ***************
 "ix_act" (area information) must be same for two conditions compared. 
  "bayes_parm.twin_meg" must be same too.

 2005/07/26 O.Y
 2005/08/03 O.Y
 2005/08/24 O.Y


 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function vb_job_chi2test_pop1(test_parm);
0002 % Chi2 test for testing significance of two time-windows of one datasets
0003 %
0004 % function vb_job_chi2test_pop1(test_parm);
0005 %
0006 % Field of 'test_parm'
0007 % ---- Input
0008 % currfile : result of VB estimation.
0009 % tsmplpre  : pre-trigger time points (time origin is 0)
0010 %             (e.g. [] (1 sample chi2 test), [1:100]')
0011 % tsmplpost : post-trigger time points (time origin is 0)
0012 %             (e.g. [1:100]', [1:300;]')
0013 %
0014 % ---- Optional Input  (default value)
0015 % testsize (0.05) : Significance level before Bonferroni correction
0016 %                   (e.g. 0.05, 0.01 and so on)
0017 % multiple_correction (ON) : Bonferroni correction is applied if ON
0018 % analyzefile ([]): MRI analyze file. If exists, the result will be
0019 %                   displayed on MRI transvese images
0020 % version   (['ANALYZE'])   : Select 'ANALYZE' or 'SBI'
0021 %                   if your vbemge version is ver0.2* or later, choose 'ANALYZE',
0022 %                   otherwise 'SBI'
0023 % outputfile ([])  : If specified, the variables 'pval','ix_dip','test_parm' are saved.
0024 %
0025 %
0026 % --- Example
0027 % test_parm.currfile = ['c00015a-a_UR4.curr.mat'];
0028 % test_parm.tsmplpre = [1:300]'; %pre-trigger
0029 % test_parm.tsmplpost = [360:420]'; %post-trigger
0030 % test_parm.analyzefile = '/home/dcfs2/oyamashi/matlab/vbmeg2/realdata/retino/SH/data/3D.img';
0031 % test_parm.version = 'ANALYZE'
0032 %
0033 % ********** Notice ***************
0034 % "ix_act" (area information) must be same for two conditions compared.
0035 %  "bayes_parm.twin_meg" must be same too.
0036 %
0037 % 2005/07/26 O.Y
0038 % 2005/08/03 O.Y
0039 % 2005/08/24 O.Y
0040 %
0041 %
0042 % Copyright (C) 2011, ATR All Rights Reserved.
0043 % License : New BSD License(see VBMEG_LICENSE.txt)
0044 
0045 vb_struct2vars(test_parm, {'currfile','tsmplpre','tsmplpost',...
0046     'testsize', 'multiple_correction',... 
0047     'analyzefile', 'version', 'outputfile'});
0048 
0049 
0050 %% default value
0051 if isempty(testsize)
0052     testsize = 0.05;
0053 end
0054 if isempty(multiple_correction)
0055     multiple_correction = ON;
0056 end
0057 
0058 %% error check
0059 if ~isempty('analyzefile')
0060     if isempty('version')
0061         fprintf('ANALYZE center coordinate system is used');
0062     end
0063 end
0064 
0065 %
0066 % SET PARAMETERS
0067 %
0068 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0069 
0070 load(currfile,'bayes_parm');
0071 % bayes_parm=change_bayes_parm_file_dir(bayes_parm,'data/');
0072 % save(bayesfile,'bayes_parm', '-append');
0073 
0074 Ntall = bayes_parm.twin_meg(2)-bayes_parm.twin_meg(1)+1;
0075 Nfile = 1;
0076 
0077 if isempty(tsmplpre)
0078     % one sample chi2 test
0079     Npop = 1;
0080     tsmpl{1} = tsmplpost;
0081 else
0082     % two sample chi2 test
0083     Npop = 2;
0084     tsmpl{1} = tsmplpost;
0085     tsmpl{2} = tsmplpre;
0086 end
0087 
0088     
0089 %
0090 % LOADING CURRENT FILES
0091 %
0092 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0093 
0094 load(currfile, 'Jact','ix_act','Varact','vb_parm');
0095 L = vb_parm.Norient;
0096 
0097 if exist('Varact') == 0
0098     error('Variance can not be found !!!');
0099 end 
0100     
0101 for i = 1 : Npop
0102     [J{i},Var{i}]=set_JV(Jact,Varact,vb_parm.Twindow,tsmpl{i});
0103 end
0104 
0105 %
0106 % TEST
0107 %
0108 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0109 [Nvact]= length(ix_act);
0110 if multiple_correction
0111     slevel = testsize/Nvact; % Bonferroni correction
0112 else
0113     slevel = testsize;
0114 end
0115     
0116 if Npop == 1
0117     pval = chi2test(J{1}, Var{1}, [], [], L);
0118     
0119     %Nt1=length(tsmpl{1});
0120     %pval = chi2test(mean(J{1},2), sum(Var{1},2)/(Nt1^2), [], [], L);
0121 else
0122     Nt1=length(tsmpl{1});
0123     Nt2=length(tsmpl{2});
0124     
0125     % same size as J{1},Var{1}
0126     Jpre = repmat(mean(J{2},2), [1,Nt1]);
0127     Vpre = repmat(sum(Var{2},2)/(Nt2^2), [1,Nt1]);
0128         
0129     pval = chi2test(J{1}, Var{1}, Jpre, Vpre, L);
0130    % pval = chi2test(mean(J{1},2), sum(Var{1},2)/(Nt1^2),...
0131    %     mean(J{2},2), sum(Var{2},2)/(Nt2^2), L);
0132 end
0133 
0134 ix_dip_tmp = find(pval > 1-slevel);
0135 ix_dip = ix_act(ix_dip_tmp);% significant vertex indices
0136 
0137 %
0138 %  PLOT THE RESULT ON MRI IMAGE
0139 %
0140 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0141 
0142 if ~isempty(analyzefile)
0143     %% read MRI image
0144     [B,Vdim,Vsize]=vb_load_analyze_to_right(analyzefile);
0145     B = reshape(B,Vdim);
0146     if strcmp(version,'SBI')
0147         load(bayes_parm.brainfile,'Vox');
0148         Vana = vb_vox_to_analyze_right(Vox,Vdim,Vsize);
0149     else
0150         V = vb_load_cortex(bayes_parm.brainfile);
0151         Vana=vb_spm_right_to_analyze_right(V,Vdim,Vsize);
0152     end
0153     vb_plot_slice(B, Vana(ix_dip,:), [90:4:190], 1, [4,4], 10);    
0154 end
0155  
0156 %
0157 % SAVE RESULT
0158 %
0159 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0160 
0161 if ~isempty(outputfile)
0162     vb_fsave(outputfile, 'ix_dip', 'pval', 'test_parm');
0163 end
0164 
0165     
0166 
0167 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0168 %
0169 % inner function
0170 %
0171 %
0172 function [J,V] = set_JV(Jact,Varact,Twindow,Tsmpl)
0173 
0174 J = Jact(:,Tsmpl);
0175 % Var for all t
0176 tstart = Tsmpl(1);
0177 tend = Tsmpl(end);
0178 for t = tstart:tend 
0179     ix_tw=find( t >= Twindow(:,1) &...
0180                 t <= Twindow(:,2));
0181     ind = t -tstart+1;
0182     Vart(:,ind) = mean(Varact(:,ix_tw),2);
0183 end
0184 V = Vart;

Generated on Tue 27-Aug-2013 11:46:04 by m2html © 2005