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)
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;