# two_group_test_spectrum

## PURPOSE

function [dz,vdz,Adz]=two_group_test_spectrum(J1,J2,p,plt,f)

## SYNOPSIS

function [dz,vdz,Adz]=two_group_test_spectrum(J1,J2,p,plt,f)

## DESCRIPTION

     function [dz,vdz,Adz]=two_group_test_spectrum(J1,J2,p,plt,f)
     Test the null hypothesis (H0) that data sets J1, J2 in
     two conditions c1,c2 have equal population spectrum

     Usage:
     [dz,vdz,Adz]=two_sample_test_spectrum(J1,J2,p)

     Inputs:
     J1   tapered fourier transform in condition 1
     J2   tapered fourier transform in condition 2
     p      p value for test (default: 0.05)
     plt    'y' for plot and 'n' for no plot
     f      frequencies (useful for plotting)

     Dimensions: J1: frequencies x number of samples in condition 1
                 J2: frequencies x number of samples in condition 2
                  number of samples = number of trials x number of tapers
     Outputs:
     dz    test statistic (will be distributed as N(0,1) under H0
     vdz   Arvesen estimate of the variance of dz
     Adz   1/0 for accept/reject null hypothesis of equal population
           coherences based dz ~ N(0,1)

     Note: all outputs are functions of frequency

     References: Arvesen, Jackkknifing U-statistics, Annals of Mathematical
     Statisitics, vol 40, no. 6, pg 2076-2100 (1969)

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SOURCE CODE

    0001 function [dz,vdz,Adz]=two_group_test_spectrum(J1,J2,p,plt,f)
    0002 % function [dz,vdz,Adz]=two_group_test_spectrum(J1,J2,p,plt,f)
    0003 % Test the null hypothesis (H0) that data sets J1, J2 in
    0004 % two conditions c1,c2 have equal population spectrum
    0005 %
    0006 % Usage:
    0007 % [dz,vdz,Adz]=two_sample_test_spectrum(J1,J2,p)
    0008 %
    0009 % Inputs:
    0010 % J1   tapered fourier transform in condition 1
    0011 % J2   tapered fourier transform in condition 2
    0012 % p      p value for test (default: 0.05)
    0013 % plt    'y' for plot and 'n' for no plot
    0014 % f      frequencies (useful for plotting)
    0015 %
    0016 %
    0017 % Dimensions: J1: frequencies x number of samples in condition 1
    0018 %             J2: frequencies x number of samples in condition 2
    0019 %              number of samples = number of trials x number of tapers
    0020 % Outputs:
    0021 % dz    test statistic (will be distributed as N(0,1) under H0
    0022 % vdz   Arvesen estimate of the variance of dz
    0023 % Adz   1/0 for accept/reject null hypothesis of equal population
    0024 %       coherences based dz ~ N(0,1)
    0025 %
    0026 %
    0027 % Note: all outputs are functions of frequency
    0028 %
    0029 % References: Arvesen, Jackkknifing U-statistics, Annals of Mathematical
    0030 % Statisitics, vol 40, no. 6, pg 2076-2100 (1969)
    0031
    0032 if nargin < 2; error('Need four sets of Fourier transforms'); end;
    0033 if nargin < 4 || isempty(plt); plt='n'; end;
    0034 %
    0035 % Test for matching dimensionalities
    0036 %
    0037 m1=size(J1,2); % number of samples, condition 1
    0038 m2=size(J2,2); % number of samples, condition 2
    0039 dof1=m1; % degrees of freedom, condition 1
    0040 dof2=m2; % degrees of freedom, condition 2
    0041 if nargin < 5 || isempty(f); f=size(J1,1); end;
    0042 if nargin < 3 || isempty(p); p=0.05; end; % set the default p value
    0043
    0044 %
    0045 % Compute the individual condition spectra, coherences
    0046 %
    0047 S1=conj(J1).*J1; % spectrum, condition 1
    0048 S2=conj(J2).*J2; % spectrum, condition 2
    0049
    0050 Sm1=squeeze(mean(S1,2)); % mean spectrum, condition 1
    0051 Sm2=squeeze(mean(S2,2)); % mean spectrum, condition 2
    0052 %
    0053 % Compute the statistic dz, and the probability of observing the value dz
    0054 % given an N(0,1) distribution i.e. under the null hypothesis
    0055 %
    0056 bias1=psi(dof1)-log(dof1); bias2=psi(dof2)-log(dof2); % bias from Thomson & Chave
    0057 var1=psi(1,dof1); var2=psi(1,dof2); % variance from Thomson & Chave
    0058 z1=log(Sm1)-bias1; % Bias-corrected Fisher z, condition 1
    0059 z2=log(Sm2)-bias2; % Bias-corrected Fisher z, condition 2
    0060 dz=(z1-z2)/sqrt(var1+var2); % z statistic
    0061
    0062 % Bug fix
    0063 % pdz=normpdf(dz,0,1); % probability of observing value dz
    0064 pdz = 2*normcdf(-abs(dz),0,1); % probability of observing value dz
    0065
    0066 %
    0067 % The remaining portion of the program computes Jackknife estimates of the mean (mdz) and variance (vdz) of dz
    0068 %
    0069 samples1=[1:m1];
    0070 samples2=[1:m2];
    0071 %
    0072 % Leave one out of one sample
    0073 %
    0074 bias11=psi(dof1-1)-log(dof1-1); var11=psi(1,dof1-1);
    0075 for i=1:m1;
    0076     ikeep=setdiff(samples1,i); % all samples except i
    0077     Sm1=squeeze(mean(S1(:,ikeep),2)); % 1 drop mean spectrum, data 1, condition 1
    0078     z1i(:,i)=log(Sm1)-bias11; % 1 drop, bias-corrected Fisher z, condition 1
    0079     dz1i(:,i)=(z1i(:,i)-z2)/sqrt(var11+var2); % 1 drop, z statistic, condition 1
    0080     ps1(:,i)=m1*dz-(m1-1)*dz1i(:,i);
    0081 end;
    0082 ps1m=mean(ps1,2);
    0083 bias21=psi(dof2-1)-log(dof2-1); var21=psi(1,dof2-1);
    0084 for j=1:m2;
    0085     jkeep=setdiff(samples2,j); % all samples except j
    0086     Sm2=squeeze(mean(S2(:,jkeep),2)); % 1 drop mean spectrum, data 2, condition 2
    0087     z2j(:,j)=log(Sm2)-bias21; % 1 drop, bias-corrected Fisher z, condition 2
    0088     dz2j(:,j)=(z1-z2j(:,j))/sqrt(var1+var21); % 1 drop, z statistic, condition 2
    0089     ps2(:,j)=m2*dz-(m2-1)*dz2j(:,j);
    0090 end;
    0091 %
    0092 % Leave one out, both samples
    0093 % and pseudo values
    0094 % for i=1:m1;
    0095 %     for j=1:m2;
    0096 %         dzij(:,i,j)=(z1i(:,i)-z2j(:,j))/sqrt(var11+var21);
    0097 %         dzpseudoval(:,i,j)=m1*m2*dz-(m1-1)*m2*dz1i(:,i)-m1*(m2-1)*dz2j(:,j)+(m1-1)*(m2-1)*dzij(:,i,j);
    0098 %     end;
    0099 % end;
    0100 %
    0101 % Jackknife mean and variance
    0102 %
    0103 % dzah=sum(sum(dzpseudoval,3),2)/(m1*m2);
    0104 ps2m=mean(ps2,2);
    0105 % dzar=(sum(ps1,2)+sum(ps2,2))/(m1+m2);
    0106 vdz=sum((ps1-ps1m(:,ones(1,m1))).*(ps1-ps1m(:,ones(1,m1))),2)/(m1*(m1-1))+sum((ps2-ps2m(:,ones(1,m2))).*(ps2-ps2m(:,ones(1,m2))),2)/(m2*(m2-1));
    0107 % vdzah=sum(sum((dzpseudoval-dzah(:,ones(1,m1),ones(1,m2))).*(dzpseudoval-dzah(:,ones(1,m1),ones(1,m2))),3),2)/(m1*m2);
    0108 %
    0109 % Test whether H0 is accepted at the specified p value
    0110 %
    0111 Adz=zeros(size(dz));
    0112 x=norminv([p/2 1-p/2],0,1);
    0113 indx=find(dz>=x(1) & dz<=x(2));
    0114 Adz(indx)=1;
    0115
    0116 % Adzar=zeros(size(dzar));
    0117 % indx=find(dzar>=x(1) & dzar<=x(2));
    0118 % Adzar(indx)=1;
    0119 %
    0120 % Adzah=zeros(size(dzah));
    0121 % indx=find(dzah>=x(1) & dzah<=x(2));
    0122 % Adzah(indx)=1;
    0123 if strcmp(plt,'y');
    0124     if isempty(f) || nargin < 5;
    0125         f=linspace(0,1,length(dz));
    0126     end;
    0127     %
    0128     % Compute the coherences
    0129     %
    0130     S1=squeeze(mean(conj(J1).*J1,2));
    0131     S2=squeeze(mean(conj(J2).*J2,2));
    0132     %
    0133     % Plot the coherence
    0134     %
    0135     subplot(311);
    0136     plot(f,S1,f,S2); legend('Data 1','Data 2');
    0137     set(gca,'FontName','Times New Roman','Fontsize', 16);
    0138     ylabel('Spectra');
    0139     title('Two group test for spectrum');
    0140     subplot(312);
    0141     plot(f,dz);
    0142     set(gca,'FontName','Times New Roman','Fontsize', 16);
    0143     ylabel('Test statistic');
    0144     conf=norminv(1-p/2,0,1);
    0145     line(get(gca,'xlim'),[conf conf]);
    0146     line(get(gca,'xlim'),[-conf -conf]);
    0147     subplot(313);
    0148     plot(f,vdz);
    0149     set(gca,'FontName','Times New Roman','Fontsize', 16);
    0150     xlabel('frequency'); ylabel('Jackknifed variance');
    0151 end;