Commit 966f2444 authored by Enrico Glerean's avatar Enrico Glerean

adding timo's version

parent 7a68d9fd
load demo_data.mat
%% ANCOVA
%Number of subjects per group
Nsubs = size(g12c1, 1);
g1_subs = 1:15;
g2_subs = 16:Nsubs;
%Extract the lower triangular of the original matrices in vector form
[~, ~, X] = find(tril(g12c1, -1));
[~, ~, Y] = find(tril(g12c2, -1));
%Regress X out of Y (the ANCOVA step). The leftover variation is in R.
[B,BINT,R] = regress(Y,[X ones(size(X))]);
%Reshape R into a lower triangular matrix in preparation for permutations
R_mat = tril(ones(Nsubs), -1);
R_mat(R_mat==1) = R;
%% Group difference via permutation test
%Extract groups 1 and 2 from R
[~, ~, R_mat_g1] = find(tril(R_mat(g1_subs, g1_subs), -1));
[~, ~, R_mat_g2] = find(tril(R_mat(g2_subs, g2_subs), -1));
%Record the t-stat value for the unpermuted voxel
[H, P, CI, STATS] = ttest2(R_mat_g1, R_mat_g2);
tstat_unpermuted = STATS.tstat;
Niter = 5000;
%initialize matrix to hold the permutations and the permuted t-stats
perms = zeros(5000, size(R_mat, 1));
tstats = zeros(Niter, 1);
%Generate random permutations
for N = 1:Niter
perms(N, :) = randperm(size(R_mat, 1))';
end
parfor N = 1:Niter
%Implement permutation on R
fake_R_mat = R_mat(perms(N, :), perms(N, :));
%Extract groups 1 and 2
[~, ~, fake_R_mat_g1] = find(tril(fake_R_mat(g1_subs, g1_subs), -1));
[~, ~, fake_R_mat_g2] = find(tril(fake_R_mat(g2_subs, g2_subs), -1));
%Calculate and record t-stat for each permutation
[H, P, CI, STATS] = ttest2(fake_R_mat_g1, fake_R_mat_g2);
tstats(N, :) = STATS.tstat;
end
%calculate p-value for unpermuted t-stat
lower_percentile = size(tstats(tstats < tstat_unpermuted), 1) / size(tstats, 1);
p_val = min(lower_percentile, 1-lower_percentile);
hist(tstats);
ylabel('Frequency');
xlabel('t-score');
hold on;
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment