Commit 8ec74acd by Enrico Glerean

adding

parent 966f2444
load demo_data.mat
%% ANCOVA
%Number of subjects per group
Nsubs = size(g12c1, 1);
g1_subs = 1:15;
g2_subs = 16:Nsubs;
G1=zeros(Nsubs);
G1(g1_subs,g1_subs)=1;
G2=zeros(Nsubs);
G2(g2_subs,g2_subs)=1;
%Extract the lower triangular of the original matrices in vector form
ids_LT=find(tril(ones(Nsubs),-1));
X = atanh(g12c1(ids_LT));
Y = atanh(g12c2(ids_LT));
%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 = ones(Nsubs);
R_mat(ids_LT) = R;
R_mat=R_mat+R_mat'; % fill in also the top triangle, used for permutations
%% Group difference via permutation test
%Extract groups 1 and 2 from R
ids_g1_LT=find(tril(G1,-1));
ids_g2_LT=find(tril(G2,-1));
R_mat_g1 = R_mat(ids_g1_LT);
R_mat_g2 = R_mat(ids_g2_LT);
%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
rng(0)
for N = 1:Niter
perms(N, :) = randperm(size(R_mat, 1))';
end
for 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=fake_R_mat(ids_g1_LT);
fake_R_mat_g2=fake_R_mat(ids_g2_LT);
%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);
[fi xi]=ksdensity(tstats,'function','cdf','npoints',200);
pval_left=interp1([-1 xi 1],[0 fi 1],tstat_unpermuted); % trick to avoid NaNs
pval=1-pval_left;
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