diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..bdb0cab --- /dev/null +++ b/.gitattributes @@ -0,0 +1,17 @@ +# Auto detect text files and perform LF normalization +* text=auto + +# Custom for Visual Studio +*.cs diff=csharp + +# Standard to msysgit +*.doc diff=astextplain +*.DOC diff=astextplain +*.docx diff=astextplain +*.DOCX diff=astextplain +*.dot diff=astextplain +*.DOT diff=astextplain +*.pdf diff=astextplain +*.PDF diff=astextplain +*.rtf diff=astextplain +*.RTF diff=astextplain diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..cd2946a --- /dev/null +++ b/.gitignore @@ -0,0 +1,47 @@ +# Windows image file caches +Thumbs.db +ehthumbs.db + +# Folder config file +Desktop.ini + +# Recycle Bin used on file shares +$RECYCLE.BIN/ + +# Windows Installer files +*.cab +*.msi +*.msm +*.msp + +# Windows shortcuts +*.lnk + +# ========================= +# Operating System Files +# ========================= + +# OSX +# ========================= + +.DS_Store +.AppleDouble +.LSOverride + +# Thumbnails +._* + +# Files that might appear in the root of a volume +.DocumentRevisions-V100 +.fseventsd +.Spotlight-V100 +.TemporaryItems +.Trashes +.VolumeIcon.icns + +# Directories potentially created on remote AFP share +.AppleDB +.AppleDesktop +Network Trash Folder +Temporary Items +.apdisk diff --git a/Example_2way_indepANOVA.m b/Example_2way_indepANOVA.m new file mode 100644 index 0000000..665e70c --- /dev/null +++ b/Example_2way_indepANOVA.m @@ -0,0 +1,94 @@ +stemFolder = pwd; +load([stemFolder filesep 'dummy'],'data'); + +cfg_neighb = []; +cfg_neighb.method = 'template'; +cfg_neighb.template = 'CTF275_neighb.mat'; +neighbours = ft_prepare_neighbours(cfg_neighb, data); + +%% Gender-Dosage test data set +% Test data set from: http://personality-project.org/r/r.guide/r.anova.html + +% Indep 2x2 ANOVA: factors Gender & dosage on alertness level + +% Design matrix for the test data +design(1,:) = [ones(1,2*4) ones(1,2*4)*2]; % First factor: Gender +design(2,:) = repmat([ones(1,4) ones(1,4)*2],1,2); % Second Factor: Dosage + +data.powspctrm = [8 12 13 12 6 7 23 14 15 12 22 14 15 12 18 22]'; + +% Results given at the personality-project.org page + +% Df Sum Sq Mean Sq F value Pr(>F) +% Gender 1 76.562 76.562 2.9518 0.1115 +% Dosage 1 5.062 5.062 0.1952 0.6665 +% Gender:Dosage 1 0.063 0.063 0.0024 0.9617 +% Residuals 12 311.250 25.938 + + + +%% Specify configuration structure: the dummy data structure consists of a single channel, time bin and frequency of interest +cfg = []; + +cfg.channel = {'MLC11'}; +cfg.latency = [1 1]; +cfg.frequency = [4 4]; +cfg.avgovertime = 'no'; +cfg.avgoverfreq = 'no'; +cfg.avgoverchan = 'no'; + +cfg.method = 'montecarlo'; +cfg.correctm = 'cluster'; +cfg.clusteralpha = 1; % consider all clusters irrespectively of significance +cfg.clusterstatistic = 'maxsum'; +cfg.minnbchan = 0; % we are looking at one channel only + +cfg.statistic = 'indepAnova2way'; +cfg.fac = 'iaxb'; % main effect Gender: 'a', main effect Dosage: 'b' and interaction: 'iaxb' +cfg.tail = 1; % F-values can only take positive values +cfg.clustertail = 1; +cfg.design = design; +cfg.neighbours = neighbours; +cfg.numrandomization = 999; % number of permutations + + + + %% define permutation strategies + main_exact_flag = 'yes'; + switch cfg.fac + case 'a' + switch main_exact_flag + case 'no' % unrestricted permutations across levels of both factors + cfg.ivar = [1 2]; + cfg.wvar = []; + cfg.uvar = []; + cfg.cvar = []; + case 'yes' % restricted permutations within levels of the other factor (exact test) + cfg.ivar = [1 2]; + cfg.wvar = []; + cfg.uvar = []; + cfg.cvar = 2; + end + + case 'b' + switch main_exact_flag + case 'no' + cfg.ivar = [1 2]; + cfg.wvar = []; + cfg.uvar = []; + cfg.cvar = []; + case 'yes' + cfg.ivar = [1 2]; + cfg.wvar = []; + cfg.uvar = []; + cfg.cvar = 1; + end + + case 'iaxb' % unrestricted permutations across levels of both factors + cfg.ivar = [1 2]; + cfg.wvar = []; + cfg.uvar = []; + cfg.cvar = []; + end + +stat = ft_freqstatistics(cfg, data); diff --git a/FtSimLink_indepANOVA.m b/FtSimLink_indepANOVA.m new file mode 100644 index 0000000..3554d88 --- /dev/null +++ b/FtSimLink_indepANOVA.m @@ -0,0 +1,69 @@ +function stat = FtSimLink_indepANOVA(data,neighbours,y,design,statfun,fac,main_exact_flag) + +data.powspctrm = y'; % data based on a dummy structure + +% Specify configuration structure +cfg = []; + +cfg.channel = {'MLC11'}; +cfg.latency = [1 1]; +cfg.frequency = [4 4]; +cfg.avgovertime = 'no'; +cfg.avgoverfreq = 'no'; +cfg.avgoverchan = 'no'; + +cfg.method = 'montecarlo'; +cfg.correctm = 'cluster'; +cfg.clusteralpha = 0.05; % threshold for a data entry being a cluster candidate +cfg.clusteralpha = 1; % show all clusters irrespectively of significance +cfg.clusterstatistic = 'maxsum'; +cfg.minnbchan = 0; % the dummy data has only one channel + +cfg.statistic = statfun; +cfg.tail = 1; % F-values can only take positive values +cfg.clustertail = 1; +cfg.design = design; +cfg.neighbours = neighbours; +cfg.numrandomization = 999; % number of permutations +cfg.fac = fac; + + % define permutations to be performed + switch cfg.fac + case 'a' % exact main a + switch main_exact_flag + case 'no' + cfg.ivar = [1 2]; + cfg.wvar = []; + cfg.uvar = []; + cfg.cvar = []; + case 'yes' + cfg.ivar = [1 2]; + cfg.wvar = []; + cfg.uvar = []; + cfg.cvar = 2; + end + + case 'b' % exact main b + switch main_exact_flag + case 'no' + cfg.ivar = [1 2]; + cfg.wvar = []; + cfg.uvar = []; + cfg.cvar = []; + case 'yes' + cfg.ivar = [1 2]; + cfg.wvar = []; + cfg.uvar = []; + cfg.cvar = 1; + end + + case 'iaxb' % a x b interaction + cfg.ivar = [1 2]; + cfg.wvar = []; + cfg.uvar = []; + cfg.cvar = []; + end + + +cfg.design = design'; +stat = ft_freqstatistics(cfg, data); \ No newline at end of file diff --git a/Readme.txt b/Readme.txt new file mode 100644 index 0000000..8c402e3 --- /dev/null +++ b/Readme.txt @@ -0,0 +1,52 @@ +How to set up Fieldtrip + +To perform a 2-way balanced, independent (between units of observation) permutation ANOVA: + +1.) put the ft_statfun_indepAnova2way.m file into the statfun folder of your current Fieldtrip version + +2.) the actual ANOVA calculations are done by the core function anova_cell_mod (adapted from Arnaud Delorme’s resampling statistical toolkit, https://de.mathworks.com/matlabcentral/fileexchange/27960-resampling-statistical-toolkit), which needs to go into the privat subfolder of the statfun folder + +3.) Next, you’ll have to modify some Fieldtrip functions. Note that the line numbers given below might be slightly different in your Fieldtrip version. They should give you an idea where to look, though. + +In ft_statistics_montecarlo: +1. around line 120 (where the defaults for the main function are set), add the following line to set the default effect of interest to be the interaction: +cfg.fac = ft_getopt(cfg, 'fac','iaxb'); + +2. starting at line 213, replace the assignment + +resample = resampledesign(cfg, design); + +with + +if isfield(cfg,'resample') + resample = cfg.resample; +else + resample = resampledesign(cfg, design); +end + +This will allow you to use some pre-computed, more complex permutation matrices (not necessary for an independent 2-way ANOVA, but e.g. for group x condition interactions in a mixed design) + +3. at approximately line 229 add + +tmpcfg.fac = cfg.fac; + +This configuration struct field will be used to indicate the factor or interaction of interest. This can be 'a' (first factor), 'b' (second factor) or 'iaxb' for the interaction + +4. then go to resampledesign (in the private folder of your Fieldtrip version) and change line 130 from + +resample = cat(2, blockres{:}); +to +resample(:,cat(2, blocksel{:})) = cat(2, blockres{:}); + +See the following bug note: http://bugzilla.fcdonders.nl/show_bug.cgi?id=1546 + +Now your Fieldtrip version is set up and ready to run permutation ANOVAs. + +You can take a look at the example script and the Sim_indepANOVA_* files (using FtSimLink_indepANOVA; you might want to decrease the number of Monte Carlo simulation to 10 in order to avoid very long running times). Before running the simulation scripts make sure that a copy of anova_cell_mod and the dummy structure are in your working directory (this won’t be necessary if you later run your own analysis in Fieldtrip). The set-up of the design matrix and of the independent and control variables cfg.ivar and cfg.cvar follows the same logic as for one-factorial designs in Fieldtrip (check out the tutorial http://www.fieldtriptoolbox.org/tutorial/cluster_permutation_timelock and http://www.fieldtriptoolbox.org/faq/how_can_i_use_the_ivar_uvar_wvar_and_cvar_options_to_precisely_control_the_permutations?s[]=design&s[]=matrix). + + +The Sim_indepANOVA_* files exemplify how to perform exact permutation tests with restricted permutations and approximate tests on either the raw data or the residuals under a reduced model as described in ‘Permutation tests for multi-factorial analysis of variance. Marti Anderson and Cajo Ter Braak. Journal of Statistical Computation and Simulation Vol. 73, Iss. 2, 2003. + +You'll have to run the permutation ANOVA for each effect separately by specifying cfg.fac = 'a', 'b' or 'iaxb'. + + diff --git a/Sim_indepANOVA_a.m b/Sim_indepANOVA_a.m new file mode 100644 index 0000000..0bbbcf5 --- /dev/null +++ b/Sim_indepANOVA_a.m @@ -0,0 +1,184 @@ +%% Simulations of an increasing main effect in an independent, balanced 2-way ANOVA without interaction +% Uses FtSimLink_indepANOVA, anova2_cell_mod + +stemFolder = pwd; +resDir = [stemFolder filesep 'SimData' filesep 'ResultsIndepANOVA_a' filesep]; +savePrefix = '2way_indep_a'; + +if ~exist(resDir,'dir') + mkdir(resDir); +end; + +Rep = 1000; % number of Monte Carlo simulations + +a = 3; % number of levels in factor A +b = 4; % number of levels in factor B +n = 5; % number of repetitions in each factor level combination + +B = [50 -50 20 -20]; +B = reshape(repmat(B,a,n),a,b,n); + +design = [repmat([1,2,3],1,b*n);repmat([1 1 1 2 2 2 3 3 3 4 4 4],1,n)]; % design matrix + + +load([stemFolder filesep 'dummy'],'data'); +cfg_neighb.method = 'template'; +cfg_neighb.template = 'CTF275_neighb.mat'; +neighbours = ft_prepare_neighbours(cfg_neighb, data); + +Method_List = {'ftest','raw','exact','res'}; +Error_List = {'exp','gauss'}; +for ee = 1:length(Error_List) + figure + hold on + for mm = 1:length(Method_List) + + method = Method_List{mm}; + err_dist = Error_List{ee}; + + if strcmp(err_dist,'exp') + sc = [10000, 40, 20, 12, 10, 8, 6, 4, 3]; + elseif strcmp(err_dist,'gauss') + sc = [100000, 200, 100, 70, 50, 40, 35, 30, 25]*3; + end + + saveStr = [savePrefix,'_',err_dist,'_',method]; + + for r = 1:length(sc) + % gradually increase effecr size of factor A + if r ==1 + A = [0; 0; 0]; + else + A = [50; 0;-50]/sc(r); + end + m_A = mean(A); + par = sqrt(sum(((A - m_A).^2)./a)); % effect size + A = reshape(repmat(A,b,n),a,b,n); + + res_A = zeros(Rep,1); + for j = 1:Rep + + % put the model together & add errors + if strcmp(err_dist,'exp') + y = A + B + (exprnd(1,a,b,n)).^3; + elseif strcmp(err_dist,'gauss') + y = A + B + randn(a,b,n); + end + + switch method + case 'ftest' % standard F-test + + % prepare data input for anova2_cell_mod + c{1,1} = squeeze(y(1,1,:))'; + c{1,2} = squeeze(y(1,2,:))'; + c{1,3} = squeeze(y(1,3,:))'; + c{1,4} = squeeze(y(1,4,:))'; + c{2,1} = squeeze(y(2,1,:))'; + c{2,2} = squeeze(y(2,2,:))'; + c{2,3} = squeeze(y(2,3,:))'; + c{2,4} = squeeze(y(2,4,:))'; + c{3,1} = squeeze(y(3,1,:))'; + c{3,2} = squeeze(y(3,2,:))'; + c{3,3} = squeeze(y(3,3,:))'; + c{3,4} = squeeze(y(3,4,:))'; + + [FA, FB, FI, dfa, dfb, dfi] = anova2_cell_mod(c); % adapted from the resampling statistic toolbox + + res_A(j) = 1 - fcdf(FA, dfa(1), dfa(2)); + + case 'raw' % permutation of raw data + fac = 'a'; + c = reshape(y,1,a*b*n); + exact = 'no'; + statfun = 'indepAnova2way'; + % the permutation ANOVA is called here + stat = FtSimLink_indepANOVA(data,neighbours,c,design,statfun,fac,exact); + res_A(j) = stat.prob; + + case 'exact' % restricted permutations + fac = 'a'; + c = reshape(y,1,a*b*n); + exact = 'yes'; + statfun = 'indepAnova2way'; + stat = FtSimLink_indepANOVA(data,neighbours,c,design,statfun,fac,exact); + res_A(j) = stat.prob; + + case 'res' % permutation of residuals + fac = 'a'; + c = reshape(y,1,a*b*n); + ncond_a = a; + ncond_b = b; + + for nfac_a = 1:ncond_a + for nfac_b = 1:ncond_b + idx_ab = design(1,:) == nfac_a & design(2,:) == nfac_b; + anovaIn{nfac_a,nfac_b} = c(idx_ab); + end + end + + + for jj = 1:size(anovaIn,2) + tmp = zeros(size(anovaIn{1,1})); + for ii = 1:size(anovaIn,1) + tmp(:,:,ii) = anovaIn{ii,jj}; + end + bmean(jj) = squeeze(mean(mean(tmp))); + end + + + for ii = 1:size(anovaIn,1) + for jj = 1:size(anovaIn,2) + idx_ab = design(1,:) == ii & design(2,:) == jj; + c_new(idx_ab) = c(idx_ab) - bmean(jj); + end + end + exact = 'no'; + statfun = 'indepAnova2way'; + stat = FtSimLink_indepANOVA(data,neighbours,c_new,design,statfun,fac,exact); + res_A(j) = stat.prob; + + end + + end + p_val(r) = length(find(res_A <= 0.05))/Rep; + paramA(r) = par; + end + + switch method + case 'ftest' + plot(paramA,p_val,':d','color',[0.4 0.4 0.4],'linewidth',2) + case 'raw' + plot(paramA,p_val,'g:d','linewidth',2) + case 'exact' + plot(paramA,p_val,':*','linewidth',2) + case 'res' + plot(paramA,p_val,'m:s','linewidth',2) + end + set(gcf,'color','w') + xlabel('\theta_A','FontSize',10,'FontName','Helvetica') + ylabel('Power','FontSize',10,'FontName','Helvetica') + if strcmp(err_dist,'gauss') + title('Gauss. errors','FontSize',11,'FontName','Helvetica','FontWeight','bold','FontAngle','italic') + elseif strcmp(err_dist,'exp') + title('Exp. errors','FontSize',11,'FontName','Helvetica','FontWeight','bold','FontAngle','italic') + end + set(gca, ... + 'Box' , 'off' , ... + 'TickDir' , 'out' , ... + 'TickLength' , [.02 .02] , ... + 'XMinorTick' , 'off' , ... + 'YMinorTick' , 'off' , ... + 'YGrid' , 'off' , ... + 'YTick' , [0.05 0.2:0.2:1], ... %'YTickLabel' , [], ... + 'LineWidth' , 1 ); + box off + xlim([min(paramA) max(paramA)]) + + cd(resDir) + save(saveStr,'paramA','p_val*','A','B','sc','method') + clear paramA p_val* c + cd(stemFolder) + end +end + + diff --git a/Sim_indepANOVA_a_spurInt.m b/Sim_indepANOVA_a_spurInt.m new file mode 100644 index 0000000..3d7559d --- /dev/null +++ b/Sim_indepANOVA_a_spurInt.m @@ -0,0 +1,191 @@ +%% Test for spurious interactions in an independent, balanced 2-way ANOVA without interaction with increasing main effect A +% Uses FtSimLink_indepANOVA, anova2_cell_mod + +stemFolder = pwd; +resDir = [stemFolder filesep 'SimData' filesep 'ResultsIndepANOVA_spurInt' filesep]; +savePrefix = '2way_indep_spurInt'; + +if ~exist(resDir,'dir') + mkdir(resDir); +end; + +Rep = 1000; % number of Monte Carlo simulations + +a = 3; % number of levels in factor A +b = 4; % number of levels in factor B +n = 5; % number of repetitions in each factor level combination + +B = [50 -50 20 -20]; +B = reshape(repmat(B,a,n),a,b,n); + +design = [repmat([1,2,3],1,b*n);repmat([1 1 1 2 2 2 3 3 3 4 4 4],1,n)]; % design matrix + +load([stemFolder filesep 'dummy'],'data'); +cfg_neighb.method = 'template'; +cfg_neighb.template = 'CTF275_neighb.mat'; +neighbours = ft_prepare_neighbours(cfg_neighb, data); + +Method_List = {'ftest','raw','res'}; +Error_List = {'exp','gauss'}; +for ee = 1:length(Error_List) + figure + hold on + for mm = 1:length(Method_List) + + method = Method_List{mm}; + err_dist = Error_List{ee}; + + % same as for Sim_indepANOVA_a + if strcmp(err_dist,'exp') + sc = [10000, 40, 20, 12, 10, 8, 6, 4, 3]; + elseif strcmp(err_dist,'gauss') + sc = [100000, 200, 100, 70, 50, 40, 35, 30, 25]*3; + end + + saveStr = [savePrefix,'_',err_dist,'_',method]; + + for r = 1:length(sc) + % gradually increase effecr size of factor A + if r ==1 + A = [0; 0; 0]; + else + A = [50; 0; -50]/sc(r); + end + m_A = mean(A); + par = sqrt(sum(((A - m_A).^2)./a)); + A = reshape(repmat(A,b,n),a,b,n); + + for j = 1:Rep + + % put the model together & add errors + if strcmp(err_dist,'exp') + y = A + B + (exprnd(1,a,b,n)).^3; + elseif strcmp(err_dist,'gauss') + y = A + B + randn(a,b,n); + end + + switch method + case 'ftest' % standard F-test + + % prepare data input for anova2_cell_mod + c{1,1} = squeeze(y(1,1,:))'; + c{1,2} = squeeze(y(1,2,:))'; + c{1,3} = squeeze(y(1,3,:))'; + c{1,4} = squeeze(y(1,4,:))'; + c{2,1} = squeeze(y(2,1,:))'; + c{2,2} = squeeze(y(2,2,:))'; + c{2,3} = squeeze(y(2,3,:))'; + c{2,4} = squeeze(y(2,4,:))'; + c{3,1} = squeeze(y(3,1,:))'; + c{3,2} = squeeze(y(3,2,:))'; + c{3,3} = squeeze(y(3,3,:))'; + c{3,4} = squeeze(y(3,4,:))'; + + + [FA, FB, FI, dfa, dfb, dfi] = anova2_cell_mod(c); % adapted from the resampling statistic toolbox + res_I(j) = 1 - fcdf(FI, dfi(1), dfi(2)); + + case 'raw' % permutation of raw data + fac = 'iaxb'; + c = reshape(y,1,a*b*n); + exact = 'no'; + statfun = 'indepAnova2way'; + % permutation ANOVA is called here + stat = FtSimLink_indepANOVA(data,neighbours,c,design,statfun,fac,exact); + res_I(j) = stat.prob; + + case 'res' % permutation of residuals + fac = 'iaxb'; + c = reshape(y,1,a*b*n); + ncond_a = a; + ncond_b = b; + + for nfac_a = 1:ncond_a + for nfac_b = 1:ncond_b + idx_ab = design(1,:) == nfac_a & design(2,:) == nfac_b; + anovaIn{nfac_a,nfac_b} = c(idx_ab); + end + end + + for ii = 1:size(anovaIn,1) + tmp = zeros(size(anovaIn{1,1})); + for jj = 1:size(anovaIn,2) + tmp(:,:,jj) = anovaIn{ii,jj}; + end + amean(ii) = squeeze(mean(mean(tmp))); + end + + for jj = 1:size(anovaIn,2) + tmp = zeros(size(anovaIn{1,1})); + for ii = 1:size(anovaIn,1) + tmp(:,:,ii) = [(anovaIn{ii,jj})]; + end + bmean(jj) = squeeze(mean(mean(tmp))); + end + + for ii = 1:size(anovaIn,1) + tmp = zeros(size(anovaIn{1,1})); + for jj = 1:size(anovaIn,2) + tmp(:,:,ii,jj) = anovaIn{ii,jj}; + end + + end + totmean = squeeze(mean(mean(mean(tmp)))); + + for ii = 1:size(anovaIn,1) + for jj = 1:size(anovaIn,2) + idx_ab = design(1,:) == ii & design(2,:) == jj; + c_new(idx_ab) = c(idx_ab) - bmean(jj) - amean(ii) + totmean; + + end + end + + exact = 'no'; + statfun = 'indepAnova2way'; + stat = FtSimLink_indepANOVA(data,neighbours,c_new,design,statfun,fac,exact); + res_I(j) = stat.prob; + + end + + end + p_val(r) = length(find(res_I <= 0.05))/Rep; + paramA(r) = par; + end + + switch method + case 'ftest' + plot(paramA,p_val,':d','color',[0.4 0.4 0.4],'linewidth',2) + case 'raw' + plot(paramA,p_val,'g:d','linewidth',2) + case 'res' + plot(paramA,p_val,'m:s','linewidth',2) + end + + set(gcf,'color','w') + xlabel('\theta_A','FontSize',10,'FontName','Helvetica') + ylabel('Power','FontSize',10,'FontName','Helvetica') + if strcmp(err_dist,'gauss') + title('Gauss. errors','FontSize',11,'FontName','Helvetica','FontWeight','bold','FontAngle','italic') + elseif strcmp(err_dist,'exp') + title('Exp. errors','FontSize',11,'FontName','Helvetica','FontWeight','bold','FontAngle','italic') + end + set(gca, ... + 'Box' , 'off' , ... + 'TickDir' , 'out' , ... + 'TickLength' , [.02 .02] , ... + 'XMinorTick' , 'off' , ... + 'YMinorTick' , 'off' , ... + 'YGrid' , 'off' , ... + 'YTick' , [0.05 0.2:0.2:1], ... + 'LineWidth' , 1 ); + box off + xlim([min(paramA) max(paramA)]) + + cd(resDir) + save(saveStr,'paramA','p_val*','A','B','sc','method') + clear paramA p_val* c + cd(stemFolder) + end +end + + diff --git a/Sim_indepANOVA_iaxb.m b/Sim_indepANOVA_iaxb.m new file mode 100644 index 0000000..3b5da27 --- /dev/null +++ b/Sim_indepANOVA_iaxb.m @@ -0,0 +1,196 @@ +%% Simulations of an increasing interaction effect in an independent, balanced 2-way ANOVA +% Uses FtSimLink_indepANOVA, anova2_cell_mod + +stemFolder = pwd; +resDir = [stemFolder filesep 'SimData' filesep 'ResultsIndepANOVA_iaxb' filesep]; +savePrefix = '2way_indep_iaxb'; + +if ~exist(resDir,'dir') + mkdir(resDir); +end; + +Rep = 1000; % number of Monte Carlo simulations + +a = 3; % number of levels in factor A +b = 4; % number of levels in factor B +n = 5; % number of repetitions in each factor level combination + +A = [50; 0; -50]; +B = [50 -50 20 -20]; +Int = A*B; +A = reshape(repmat(A,b,n),a,b,n); +B = reshape(repmat(B,a,n),a,b,n); + +design = [repmat([1,2,3],1,b*n);repmat([1 1 1 2 2 2 3 3 3 4 4 4],1,n)]; + +load([stemFolder filesep 'dummy'],'data'); + +cfg_neighb.method = 'template'; +cfg_neighb.template = 'CTF275_neighb.mat'; +neighbours = ft_prepare_neighbours(cfg_neighb, data); + +Method_List = {'ftest','raw','res'}; +Error_List = {'exp','gauss'}; +for ee = 1:length(Error_List) + figure + hold on + for mm = 1:length(Method_List) + + method = Method_List{mm}; + err_dist = Error_List{ee}; + + if strcmp(err_dist,'exp') + sc = [100000, 200, 100, 80, 60, 50, 40, 35, 30]*3; + elseif strcmp(err_dist,'gauss') + sc = [100000, 200, 100, 70, 50, 40, 30, 25, 20]*100; + end + saveStr = [savePrefix,'_',err_dist,'_',method]; + + for r = 1:length(sc) + if r ==1 + AB = Int/sc(r); + AB = zeros(size(AB)); + else + AB = Int/sc(r); + end + % gradually increase effect size of interaction factor + m_AB = mean(mean(AB)); + par = sqrt(sum(sum(((AB - m_AB).^2)./(a*b)))); + AB = reshape(repmat(AB,1,n),a,b,n); + + for j = 1:Rep + + % put the model together & add errors + if strcmp(err_dist,'exp') + y = A + B + AB + (exprnd(1,a,b,n)).^3; + elseif strcmp(err_dist,'gauss') + y = A + B + AB + randn(a,b,n); + end + + switch method + case 'ftest' % standard F-test + + % prepare data input for anova2_cell_mod + c{1,1} = squeeze(y(1,1,:))'; + c{1,2} = squeeze(y(1,2,:))'; + c{1,3} = squeeze(y(1,3,:))'; + c{1,4} = squeeze(y(1,4,:))'; + c{2,1} = squeeze(y(2,1,:))'; + c{2,2} = squeeze(y(2,2,:))'; + c{2,3} = squeeze(y(2,3,:))'; + c{2,4} = squeeze(y(2,4,:))'; + c{3,1} = squeeze(y(3,1,:))'; + c{3,2} = squeeze(y(3,2,:))'; + c{3,3} = squeeze(y(3,3,:))'; + c{3,4} = squeeze(y(3,4,:))'; + + + [FA, FB, FI, dfa, dfb, dfi] = anova2_cell_mod(c); % adapted from the resampling statistic toolbox + + res_I(j) = 1 - fcdf(FI, dfi(1), dfi(2)); + + case 'raw' % permutation of raw data + fac = 'iaxb'; + c = reshape(y,1,a*b*n); + exact = 'no'; + statfun = 'indepAnova2way'; + % permutation ANOVA is called here + stat = FtSimLink_indepANOVA(data,neighbours,c,design,statfun,fac,exact); + res_I(j) = stat.prob; + + + case 'res' % permutation of residuals + fac = 'iaxb'; + c = reshape(y,1,a*b*n); + ncond_a = a; + ncond_b = b; + + for nfac_a = 1:ncond_a + for nfac_b = 1:ncond_b + idx_ab = design(1,:) == nfac_a & design(2,:) == nfac_b; + anovaIn{nfac_a,nfac_b} = c(idx_ab); + end + end + + for ii = 1:size(anovaIn,1) + tmp = zeros(size(anovaIn{1,1})); + for jj = 1:size(anovaIn,2) + tmp(:,:,jj) = anovaIn{ii,jj}; + end + amean(ii) = squeeze(mean(mean(tmp))); + end + + for jj = 1:size(anovaIn,2) + tmp = zeros(size(anovaIn{1,1})); + for ii = 1:size(anovaIn,1) + tmp(:,:,ii) = [(anovaIn{ii,jj})]; + end + bmean(jj) = squeeze(mean(mean(tmp))); + end + + for ii = 1:size(anovaIn,1) + tmp = zeros(size(anovaIn{1,1})); + for jj = 1:size(anovaIn,2) + tmp(:,:,ii,jj) = anovaIn{ii,jj}; + end + + end + totmean = squeeze(mean(mean(mean(tmp)))); + + for ii = 1:size(anovaIn,1) + for jj = 1:size(anovaIn,2) + idx_ab = design(1,:) == ii & design(2,:) == jj; + c_new(idx_ab) = c(idx_ab) - bmean(jj) - amean(ii) + totmean; + + end + end + + exact = 'no'; + statfun = 'indepAnova2way'; + stat = FtSimLink_indepANOVA(data,neighbours,c_new,design,statfun,fac,exact); + res_I(j) = stat.prob; + + end + + end + p_val(r) = length(find(res_I <= 0.05))/Rep; + paramAB(r) = par; + end + + %% + switch method + case 'ftest' + plot(paramAB,p_val,':d','color',[0.4 0.4 0.4],'linewidth',2) + case 'raw' + plot(paramAB,p_val,'g:d','linewidth',2) + case 'res' + plot(paramAB,p_val,'m:s','linewidth',2) + end + + set(gcf,'color','w') + xlabel('\theta_{AB}','FontSize',10,'FontName','Helvetica') + ylabel('Power','FontSize',10,'FontName','Helvetica') + if strcmp(err_dist,'gauss') + title('Gauss. errors','FontSize',11,'FontName','Helvetica','FontWeight','bold','FontAngle','italic') + elseif strcmp(err_dist,'exp') + title('Exp. errors','FontSize',11,'FontName','Helvetica','FontWeight','bold','FontAngle','italic') + end + set(gca, ... + 'Box' , 'off' , ... + 'TickDir' , 'out' , ... + 'TickLength' , [.02 .02] , ... + 'XMinorTick' , 'off' , ... + 'YMinorTick' , 'off' , ... + 'YGrid' , 'off' , ... + 'YTick' , [0.05 0.2:0.2:1], ... + 'LineWidth' , 1 ); + box off + xlim([min(paramAB) max(paramAB)]) + cd(resDir) + save(saveStr,'paramAB','p_val*','A','B', 'AB', 'sc','method') + clear paramAB p_val* c + cd(stemFolder) + end +end + + diff --git a/anova2_cell_mod.m b/anova2_cell_mod.m new file mode 100644 index 0000000..2fcb961 --- /dev/null +++ b/anova2_cell_mod.m @@ -0,0 +1,109 @@ +%% Modified version of the ANOVA2_CELL function of Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 2010 +% Source: http://www.mathworks.com/matlabcentral/fileexchange/27960-resampling-statistical-toolkit/content/statistics/anova2_cell.m +% +% Header A. Delorme +% anova2_cell() - compute F-values in cell array using ANOVA. +% +% Usage: +% >> [FC FR FI dfc dfr dfi] = anova2_cell( data ); +% +% Inputs: +% data = data consisting of PAIRED arrays to be compared. The last +% dimension of the data array is used to compute ANOVA. +% Outputs: +% FC - F-value for columns. +% FR - F-value for rows. +% FI - F-value for interaction. +% dfc - degree of freedom for columns. +% dfr - degree of freedom for rows. +% dfi - degree of freedom for interaction. +% +% Note: the advantage over the ANOVA2 function of Matlab statistical +% toolbox is that this function works on arrays (see examples). Note +% also that you still need the statistical toolbox to assess +% significance using the fcdf() function. The other advantage is that +% this function will work with complex numbers. +% +% Example: +% a = { rand(1,10) rand(1,10) rand(1,10); rand(1,10) rand(1,10) rand(1,10) } +% [FC FR FI dfc dfr dfi] = anova2_cell(a) +% signifC = 1-fcdf(FC, dfc(1), dfc(2)) +% signifR = 1-fcdf(FR, dfr(1), dfr(2)) +% signifI = 1-fcdf(FI, dfi(1), dfi(2)) +% +% % for comparison +% anova2( [ a{1,1}' a{1,2}' a{1,3}'; a{2,1}' a{2,2}' a{2,3}' ], 10) +% +% b = { [ a{1,1}; a{1,1} ] [ a{1,2}; a{1,2} ] [ a{1,3}; a{1,3} ]; +% [ a{2,1}; a{2,1} ] [ a{2,2}; a{2,2} ] [ a{2,3}; a{2,3} ] } +% [FC FR FI dfc dfr dfi] = anova2_cell(b) +% +% c{1,1} = reshape(repmat(b{1,1}, [2 1]),2,2,10); +% c{1,2} = reshape(repmat(b{1,2}, [2 1]),2,2,10); +% c{1,3} = reshape(repmat(b{1,3}, [2 1]),2,2,10); +% c{2,3} = reshape(repmat(b{2,3}, [2 1]),2,2,10); +% c{2,2} = reshape(repmat(b{2,2}, [2 1]),2,2,10); +% c{2,1} = reshape(repmat(b{2,1}, [2 1]),2,2,10) +% [FC FR FI dfc dfr dfi] = anova2_cell(c) +% +% Author: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 2005 +% +% Reference: +% Schaum's outlines in statistics (3rd edition). 1999. Mc Graw-Hill. + +% Copyright (C) Arnaud Delorme +% +% This program is free software; you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation; either version 2 of the License, or +% (at your option) any later version. +% +% This program is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program; if not, write to the Free Software +% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +function [FA, FB, FI, freeA, freeB, freeI] = anova2_cell_mod(data) + +% compute all means and all std +% ----------------------------- +a = size(data,1); +b = size(data,2); +c = size(data{1}, 2); + + + +VE = zeros( size(data{1},1),1); +m = zeros( [ size(data{1},1) size(data) ] ); +for i = 1:a + for ii = 1:b + tmpm = mean(data{i,ii}, 2); + m(:,i,ii) = tmpm; + VE = VE+sum( (data{i,ii}-repmat(tmpm, [1 size(data{i,ii},2)])).^2, 2); + end; +end; +X = mean(mean(m,3),2); +Xj = mean(m,3); +Xk = mean(m,2); +VR = b*c*sum( (Xj-repmat(X, [1 size(Xj,2)])).^2, 2 ); +VC = a*c*sum( (Xk-repmat(X, [1 1 size(Xk,3)])).^2, 3 ); + +Xj = repmat(Xj, [1 1 size(m,3) ]); +Xk = repmat(Xk, [1 size(m,2) 1]); +VI = c*sum( sum( ( m - Xj - Xk + repmat(X, [1 size(m,2) size(m,3)]) ).^2, 3), 2 ); +SR2 = VR/(a-1); +SC2 = VC/(b-1); +SI2 = VI/(a-1)/(b-1); +SE2 = VE/(a*b*(c-1)); + +FA = SR2./SE2; % rows +FB = SC2./SE2; % columns +FI = SI2./SE2; % interaction + +freeA = [ a-1 a*b*(c-1) ]; +freeB = [ b-1 a*b*(c-1) ]; +freeI = [ (a-1)*(b-1) a*b*(c-1) ]; \ No newline at end of file diff --git a/dummy.mat b/dummy.mat new file mode 100644 index 0000000..fec94d3 Binary files /dev/null and b/dummy.mat differ diff --git a/ft_statfun_indepAnova2way.m b/ft_statfun_indepAnova2way.m new file mode 100644 index 0000000..c898d1b --- /dev/null +++ b/ft_statfun_indepAnova2way.m @@ -0,0 +1,169 @@ +function [s,cfg] = ft_statfun_indepAnova2way(cfg, dat, design) + +% FT_STATFUN_INDEPANOVA2WAY calculates the independent samples 2-way ANOVA +% F-statistic(s) +% on the biological data in dat (the dependent variable), using the information on +% the independent variables (iv) in design. +% +% Use this function by calling one of the high-level statistics functions as: +% [stat] = ft_timelockstatistics(cfg, timelock1, timelock2, ...) +% [stat] = ft_freqstatistics(cfg, freq1, freq2, ...) +% [stat] = ft_sourcestatistics(cfg, source1, source2, ...) +% with the following configuration option: +% cfg.statistic = 'indepAnova2way' +% see FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS or FT_SOURCESTATISTICS for details. +% +% For low-level use, the external interface of this function has to be +% [s,cfg] = ft_statfun_indepAnova2way(cfg, dat, design); +% where +% dat contains the biological data, Nsamples x Nreplications +% design contains the independent variable (iv), Nfac x Nreplications +% +% Configuration options: +% cfg.computestat = 'yes' or 'no', calculate the statistic (default='yes') +% cfg.computecritval = 'yes' or 'no', calculate the critical values of the test statistics (default='no') +% cfg.computeprob = 'yes' or 'no', calculate the p-values (default='no') +% +% The following options are relevant if cfg.computecritval='yes' and/or +% cfg.computeprob='yes'. +% cfg.alpha = critical alpha-level of the statistical test (default=0.05) +% cfg.tail = -1, 0, or 1, left, two-sided, or right (default=1) +% cfg.tail in combination with cfg.computecritval='yes' +% determines whether the critical value is computed at +% quantile cfg.alpha (with cfg.tail=-1), at quantiles +% cfg.alpha/2 and (1-cfg.alpha/2) (with cfg.tail=0), or at +% quantile (1-cfg.alpha) (with cfg.tail=1). +% +% Design specification: +% cfg.ivar = row number(s) of the design that contain(s) the labels of the conditions that must be +% compared (default= [1,2]). The labels range from 1 to the number of conditions. + +% Adapted by Saskia Helbling from Eric Maris' statfun_depsamplesT +% Independent two-way ANOVA is calculated by means of the function +% ANOVA2_CELL_MOD based on ANOVA2_CELL from the Resampling statistical toolkit by Arnaud Delorme +% http://www.mathworks.com/matlabcentral/fileexchange/27960-resampling-statistical-toolkit/content/statistics/anova2_cell.m +% published under the GNU General Public License as published by +% the Free Software Foundation; either version 2 of the License, or +% (at your option) any later version. + +% Subversion does not use the Log keyword, use 'svn log ' or 'svn -v log | less' to get detailled information + +% set the defaults +if ~isfield(cfg, 'computestat'), cfg.computestat='yes'; end; +if ~isfield(cfg, 'computecritval'), cfg.computecritval='no'; end; +if ~isfield(cfg, 'computeprob'), cfg.computeprob='no'; end; +if ~isfield(cfg, 'alpha'), cfg.alpha=0.05; end; +if ~isfield(cfg, 'tail'), cfg.tail=1; end; + +%% perform some checks on the configuration +if strcmp(cfg.computeprob,'yes') && strcmp(cfg.computestat,'no') + error('P-values can only be calculated if the test statistics are calculated.'); +end; +if isfield(cfg,'uvar') && ~isempty(cfg.uvar) + error('cfg.uvar should not exist for an independent samples statistic'); +end + +% only balanced 2-way factorial ANOVA for now +ncond_a = length(unique(design(cfg.ivar(1),:))); +ncond_b = length(unique(design(cfg.ivar(2),:))); + +ncond = ncond_a*ncond_b; +nrepl=zeros(ncond_a,ncond_b); + +for condindx_a=1:ncond_a + for condindx_b=1:ncond_b + nrepl(condindx_a,condindx_b)=nrepl(condindx_a,condindx_b)+length(find(design(cfg.ivar(1),:)==condindx_a& design(cfg.ivar(2),:)==condindx_b)); + end +end; +if min(nrepl)~=max(nrepl) + error('Use only balanced designs for 2-way (or higher order) ANOVA!'); % this condition might be relaxed n the future +end + +if sum(sum(nrepl))