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