Nsim = 10000;%%% number of simulations
Nsub = 20;%%% number of participants in each group
Nvar = 20;%%% number of variables that you measure
%% two repliactions, only Nsub per group
for k=1:Nsim,
% random variables - two studies A and B
A= randn(Nsub,Nvar);
B= randn(Nsub,Nvar);
% correlation between each variable for each of the two studies
% separately
[Ra,pa]=corrcoef(A);
[Rb,pb]=corrcoef(B);
% extracting the relevant p-values and correlations (upper triangular matrix without diagonal)
Ta = triu(pa,1);
Tb = triu(pb,1);
Rpa = triu(Ra,1);
Rpb = triu(Rb,1);
% matrix to vectors for non-zero entries
Ca = Ta(Ta~=0);
Cb = Tb(Tb~=0);
Rca = Rpa(Ta~=0);
Rcb = Rpb(Tb~=0);
% detecting significant correlations, no corrections for multiple comparisons
PVa = Ca<0.05;
PVb = Cb<0.05;
% detecting significant correlations, with Bonferroni corrections for multiple comparisons
Ncorr = (Nvar^2 - Nvar)/2;
PVawc = Ca<0.05/Ncorr;
PVbwc = Cb<0.05/Ncorr;
% checking that the correlations have the same sign
Sab = sign(Rca).*sign(Rcb)>0;
% detecting when the two same variables correlate in the two
% replications
Sig2(k) = sum(PVa.*PVb.*Sab);
Sig2wc(k) = sum(PVawc.*PVbwc.*Sab);
end
% average number of correlated pairs dectected across the two replications
% Bear in mind, there should be none...
disp(['average number of significant correlations present in both replications: ' num2str(mean(Sig2))])
disp(['average number of significant correlations present in both replications, corrected for mult. comp.: ' num2str(mean(Sig2wc))])
%% one experiment with 2*Nsub per group
for k=1:Nsim,
% random variables
C= randn(2*Nsub,Nvar);
% correlation between each variable
[Rc,pc]=corrcoef(C);
% extracting the relevant p-values (upper triangular matrix without diagonal)
Tc = triu(pc,1);
% matrix to vectors for non-zero entries
Cc = Tc(Tc~=0);
% detecting significant correlations without mult. comp. corr
PVc = Cc<0.05;
% detecting significant correlations with mult. comp. corr
Ncorr = (Nvar^2 - Nvar)/2;
PVcwc = Cc<0.05/Ncorr;
% detecting when the two same variables correlate in the two
% replications
Sig1(k) = sum(PVc);
Sig1wc(k) = sum(PVcwc);
end
% average number significant correlations detected
% Bear in mind, there should be none...
disp(['average number of significant correlations present in one bigger studies: ' num2str(mean(Sig1))])
disp(['average number of significant correlations present in both replications, corrected for mult. comp.: ' num2str(mean(Sig1wc))])