function p=scrambledpvalue(A,B) % scrambledpvalue calculates the negative log of the p value when % correlating two input signals A and B. A is a two-dimensional signal % (e.g. gamma band power over all channels) and B is a one-dimensional % signal (e.g. speech power). The p-values are computed using a % randomization approach, where A is flipped, circularly shifted by a % random amount and correlated with B for a large number of simulations % (default being 500). Then the distribution of r-values found from % these simulations are fit to a beta distribution to find the % parameters of this random distribution. The p-value is then found by % finding where the true r-value from correlating A with B lies on this % random distribution. % CALLING SEQUENCE: % scambledpvalue(A,B) % INPUT: % A: Two-dimensional signal (e.g. gamma band across all ECoG channels) % that we want to correlate with B % B: One-dimensional signal (e.g. speech) or two-dimensional signal % (e.g. covert gamma band power) that we want to correlate with A. % Specify the number of simulations. Default is 500. numsim=100; % Flip A and circularly shift it by a random number of samples. Repeat this % for the number of simulations defined above. Populate a matrix 'data' % with the scaled r-values for each iteration. d=flipud(A); data=zeros(numsim,size(A,2)); for i=1:numsim r=randperm(size(A,1)); g=circshift(d,r(1)); parfor j=1:size(A,2) if size(B,2)==1 data(i,j)=(corr(g(:,j),B(:,1))/2)+0.5; else data(i,j)=(corr(g(:,j),B(:,j))/2)+0.5; end end end %Set the nan values to zero. ind=isnan(data); index= ind(1,:)==1; data(:,index)=0; % Fit a beta distribution to the random distribution of scaled r-values % obtained above. Find p-value using the true scaled r-value from % correlating A and B and fitting it to this random beta distribution. p=zeros(1,size(A,2)); for i=1:size(data,2) if length(unique(data(:,i)))~=1 param=mle(data(:,i),'distribution','beta'); if size(B,2)==1 p(1,i)=1-betacdf((corr(A(:,i),B(:,1))/2)+0.5,param(1),param(2)); else p(1,i)=1-betacdf((corr(A(:,i),B(:,i))/2)+0.5,param(1),param(2)); end else p(1,i)=1; end end % Add an epsilon value (10^-20) to the p-values to get rid of zeros prior % to finding the -log of the p-values. These -log p values are then % returned to the calling function/script. p=p'; p=p+(10^-20); p=-log(p); end