|
| 1 | +% Taku Ito |
| 2 | +% 6/16/2017 |
| 3 | + |
| 4 | +% Test script to test code for permutation testing in matlab |
| 5 | +clear all |
| 6 | +clc |
| 7 | + |
| 8 | +addpath('../matlabCode/') |
| 9 | +addpath('../matlabCode/mult_comp_perm_t1/') |
| 10 | + |
| 11 | +permutations=10000; |
| 12 | + |
| 13 | +%% Create data set |
| 14 | +numVoxels = 10000; |
| 15 | +numSubjs = 20; |
| 16 | +nCond = 2; |
| 17 | +sigEffects = 20; |
| 18 | +effectAmp = 3; |
| 19 | + |
| 20 | +disp(['Running simulation analysis with ' num2str(permutations) ' permutations...']) |
| 21 | +disp(['Running contrasts on ' num2str(numVoxels) ' voxels (' num2str(numVoxels) ' independent tests) with ' num2str(numSubjs) ' subjects']) |
| 22 | +disp([num2str(nCond) ' conditions, with ' num2str(sigEffects) ' significant real effets']) |
| 23 | + |
| 24 | +% Construct data set (already contrasted) |
| 25 | +dataSet = zeros(numVoxels,numSubjs,nCond); |
| 26 | +for vox=1:numVoxels |
| 27 | + if vox <= sigEffects |
| 28 | + dataSet(vox,:, 1) = randn(1,numSubjs) + effectAmp; |
| 29 | + dataSet(vox,:, 2) = randn(1,numSubjs); |
| 30 | + else |
| 31 | + dataSet(vox,:, 1) = randn(1,numSubjs); |
| 32 | + dataSet(vox,:, 2) = randn(1,numSubjs); |
| 33 | + end |
| 34 | +end |
| 35 | + |
| 36 | +% Create contrast map |
| 37 | +contrastSet = dataSet(:,:,1) - dataSet(:,:,2); |
| 38 | + |
| 39 | + |
| 40 | +%% Run permutation test |
| 41 | +disp('Timing permutation testing using House code') |
| 42 | +tic |
| 43 | +alpha = .05; |
| 44 | +[t, maxT_thresh] = maxT(contrastSet,'nullmean', 0, ... |
| 45 | + 'alpha',alpha,... |
| 46 | + 'tail', 1, ... |
| 47 | + 'permutations',permutations, ... |
| 48 | + 'nproc', 20); |
| 49 | +toc |
| 50 | + |
| 51 | +%disp(['Number of true effects: ' num2str(sigEffects)]) |
| 52 | +%disp(['Number of statistically significant effects (p < .05): ' num2str(sum(p_fwe>.95))]) |
| 53 | + |
| 54 | + |
| 55 | +%% Compare with Groppe's Matlab function |
| 56 | +% Format data in his required format |
| 57 | +tmpdat = contrastSet'; |
| 58 | +disp('Running permutations using Groppe''s Matlab function') |
| 59 | +tic |
| 60 | +[pval, t_orig, crit_t, est_alpha, seed_state]=mult_comp_perm_t1(tmpdat,permutations,1,alpha,0,0); |
| 61 | +toc |
0 commit comments