diff --git a/permutools/booteffectsize.m b/permutools/booteffectsize.m index 04f7a4e..4b88dcc 100644 --- a/permutools/booteffectsize.m +++ b/permutools/booteffectsize.m @@ -1,4 +1,4 @@ -function [d,ci,stats,bdist] = booteffectsize(x,m,varargin) +function [d,ci,stats,dist] = booteffectsize(x,m,varargin) %BOOTEFFECTSIZE Effect size with bootstrapped confidence intervals. % D = BOOTEFFECTSIZE(X) returns the effect size measure for a single % sample X based on Cohen's d. By default, Cohen's d is bias-corrected @@ -33,8 +33,8 @@ % 'sd' -- the pooled standard deviation, or of X for a one- % sample or Glass' delta measure % -% [D,CI,STATS,BDIST] = BOOTEFFECTSIZE(...) returns the bootstrap -% distribution of the effect size measure. +% [D,CI,STATS,DIST] = BOOTEFFECTSIZE(...) returns the bootstrapped +% sampling distribution of the effect size measure. % % [...] = BOOTEFFECTSIZE(...,'PARAM1',VAL1,'PARAM2',VAL2,...) specifies % additional parameters and their values. Valid parameters are the @@ -93,12 +93,14 @@ % PERMUTOOLS https://github.com/mickcrosse/PERMUTOOLS % References: -% [1] Hentschke H, Stuttgen MC (2011) Computation of measures of +% [1] Crosse MJ, Foxe JJ, Molholm S (2024) PERMUTOOLS: A MATLAB +% Package for Multivariate Permutation Testing. arXiv 2401.09401. +% [2] Hentschke H, Stuttgen MC (2011) Computation of measures of % effect size for neuroscience data sets. Eur J Neurosci, % 34:1887–1894. -% [2] Cohen J (1969) Statistical power for the behavioural sciences. +% [3] Cohen J (1969) Statistical power for the behavioural sciences. % London: Academic Press. -% [3] Hedges LV, Olkin I (1985) Statistical methods for meta- +% [4] Hedges LV, Olkin I (1985) Statistical methods for meta- % analysis. San Diego, CA: Academic Press. % © 2018-2023 Mick Crosse @@ -289,13 +291,12 @@ if nargout > 1 - % Estimate bootstrap distribution & CIs rng(arg.seed); ci = zeros(2,nvar); for i = 1:nvar - % Bootstrap samples + % Generate random bootstraps xb = x(:,i); yb = y(:,i); idx = ceil(nobsx(i)*rand([nobsx(i),arg.nboot])); @@ -305,15 +306,12 @@ end yb = yb(idx); - % Compute sample variance using fast algo + % Estimate sampling distribution smxb = sum(xb,nanflag); smyb = sum(yb,nanflag); varxb = (sum(xb.^2,nanflag)-(smxb.^2)/nobsx(i))/dfx(i); varyb = (sum(yb.^2,nanflag)-(smyb.^2)/nobsy(i))/dfy(i); - if arg.paired - - % Compute difference between samples switch arg.effect case 'cliff' diffxyb = zeros(max(nobsx)^2,arg.nboot); @@ -325,21 +323,14 @@ otherwise diffxyb = xb-yb; end - - % Compute mean difference switch arg.effect case 'mediandiff' mub = median(diffxyb,nanflag); otherwise mub = sum(diffxyb,nanflag)/nobs(i); end - - % Compute standard deviation sdb = sqrt((varxb+varyb)/2); - else - - % Compute difference between samples switch arg.effect case 'cliff' diffxyb = zeros(max(nobsx)*max(nobsy),arg.nboot); @@ -348,8 +339,6 @@ diffxyb(:,j) = diffi(:); end end - - % Compute mean difference switch arg.effect case 'cliff' mub = sum(diffxyb,nanflag)/nobs(i); @@ -358,8 +347,6 @@ otherwise mub = smxb/nobsx(i)-smyb/nobsy(i); end - - % Compute standard deviation switch arg.effect case 'glass' sdb = sqrt(varxb); @@ -371,19 +358,16 @@ sdb = sqrt((varxb+varyb)/2); end end - end - - % Compute effect size switch arg.effect case {'cliff','meandiff','mediandiff'} - bdist = mub; + dist = mub; otherwise - bdist = mub./sdb; + dist = mub./sdb; end - % Compute confidence intervals - ci(:,i) = prctile(bdist,100*[arg.alpha/2;1-arg.alpha/2]); + % Compute CI + ci(:,i) = prctile(dist,100*[arg.alpha/2;1-arg.alpha/2]); end diff --git a/permutools/permuanova1.m b/permutools/permuanova1.m index 75ab8db..9d4cbe6 100644 --- a/permutools/permuanova1.m +++ b/permutools/permuanova1.m @@ -1,4 +1,4 @@ -function [f,p,ci,stats,tbl,pdist] = permuanova1(x,group,varargin) +function [f,p,ci,stats,tbl,dist] = permuanova1(x,group,varargin) %PERMUANOVA1 One-way permutation-based analysis of variance (ANOVA). % F = PERMUANOVA1(X) performs a one-way permutation-based ANOVA for % comparing the means of two or more groups of data in matrix X, and @@ -35,8 +35,8 @@ % [F,P,CI,STATS,TBL] = PERMUANOVA1(...) returns the ANOVA table contents % as a cell array. % -% [F,P,CI,STATS,TBL,PDIST] = PERMUANOVA1(...) returns the permutation -% distribution of the test statistic. +% [F,P,CI,STATS,TBL,DIST] = PERMUANOVA1(...) returns the permuted +% sampling distribution of the test statistic. % % [...] = PERMUANOVA1(...,'PARAM1',VAL1,'PARAM2',VAL2,...) specifies % additional parameters and their values. Valid parameters are the @@ -60,13 +60,8 @@ % PERMUTOOLS https://github.com/mickcrosse/PERMUTOOLS % References: -% [1] Blair RC, Higgins JJ, Karniski W, Kromrey JD (1994) A Study of -% Multivariate Permutation Tests Which May Replace Hotelling's T2 -% Test in Prescribed Circumstances. Multivariate Behav Res, -% 29(2):141-163. -% [2] Groppe DM, Urbach TP, Kutas M (2011) Mass univariate analysis -% of event-related brain potentials/fields I: A critical tutorial -% review. Psychophysiology, 48(12):1711-1725. +% [1] Crosse MJ, Foxe JJ, Molholm S (2024) PERMUTOOLS: A MATLAB +% Package for Multivariate Permutation Testing. arXiv 2401.09401. % © 2018-2024 Mick Crosse % CNL, Albert Einstein College of Medicine, NY. @@ -140,8 +135,8 @@ maxnobs = numel(x); [~,idx] = sort(rand(maxnobs,arg.nperm)); - % Estimate permutation distribution - pdist = zeros(arg.nperm,1); + % Estimate sampling distribution + dist = zeros(arg.nperm,1); for i = 1:arg.nperm xp = x(idx(:,i)); xp = reshape(xp,shapex); @@ -153,13 +148,13 @@ essp = essp + sum((xg-mp).^2,'all',nanflag); rssp = rssp + np*sum((mp-gm).^2,'all',nanflag); end - pdist(i) = (rssp/dfr)/(essp/dfe); + dist(i) = (rssp/dfr)/(essp/dfe); end - % Compute p-value and CIs - p = (sum(f<=pdist)+1)/(arg.nperm+1); + % Compute p-value & CI + p = (sum(f<=dist)+1)/(arg.nperm+1); if nargout > 2 - crit = prctile(pdist,100*(1-arg.alpha)); + crit = prctile(dist,100*(1-arg.alpha)); ci = [f./crit;Inf]; end diff --git a/permutools/permuanova2.m b/permutools/permuanova2.m index e118c4e..4b09ef8 100644 --- a/permutools/permuanova2.m +++ b/permutools/permuanova2.m @@ -1,4 +1,4 @@ -function [f,p,ci,stats,tbl,pdist] = permuanova2(x,reps,varargin) +function [f,p,ci,stats,tbl,dist] = permuanova2(x,reps,varargin) %PERMUANOVA2 Two-way permutation-based analysis of variance (ANOVA). % F = PERMUANOVA2(X) performs a balanced two-way permutation-based ANOVA % for comparing the means of two or more columns and two or more rows of @@ -36,8 +36,8 @@ % [F,P,CI,STATS,TBL] = PERMUANOVA2(...) returns the ANOVA table contents % as a cell array. % -% [F,P,CI,STATS,TBL,PDIST] = PERMUANOVA2(...) returns the permutation -% distributions of the test statistics. +% [F,P,CI,STATS,TBL,PDIST] = PERMUANOVA2(...) returns the permuted +% sampling distributions of the test statistics. % % [...] = PERMUANOVA2(...,'PARAM1',VAL1,'PARAM2',VAL2,...) specifies % additional parameters and their values. Valid parameters are the @@ -61,13 +61,8 @@ % PERMUTOOLS https://github.com/mickcrosse/PERMUTOOLS % References: -% [1] Blair RC, Higgins JJ, Karniski W, Kromrey JD (1994) A Study of -% Multivariate Permutation Tests Which May Replace Hotelling's T2 -% Test in Prescribed Circumstances. Multivariate Behav Res, -% 29(2):141-163. -% [2] Groppe DM, Urbach TP, Kutas M (2011) Mass univariate analysis -% of event-related brain potentials/fields I: A critical tutorial -% review. Psychophysiology, 48(12):1711-1725. +% [1] Crosse MJ, Foxe JJ, Molholm S (2024) PERMUTOOLS: A MATLAB +% Package for Multivariate Permutation Testing. arXiv 2401.09401. % © 2018-2024 Mick Crosse % CNL, Albert Einstein College of Medicine, NY. @@ -177,11 +172,11 @@ maxnobs = numel(x); [~,idx] = sort(rand(maxnobs,arg.nperm)); - % Estimate permutation distributions - pdistc = zeros(arg.nperm,1); - pdistr = zeros(arg.nperm,1); + % Estimate sampling distributions + distc = zeros(arg.nperm,1); + distr = zeros(arg.nperm,1); if reps > 1 - pdisti = zeros(arg.nperm,1); + disti = zeros(arg.nperm,1); end for i = 1:arg.nperm xp = x(idx(:,i)); @@ -203,18 +198,18 @@ essp = issp; end msep = essp/dfe; - pdistc(i) = cssp/dfc/msep; - pdistr(i) = rssp/dfr/msep; + distc(i) = cssp/dfc/msep; + distr(i) = rssp/dfr/msep; if reps > 1 - pdisti(i) = issp/dfi/msep; + disti(i) = issp/dfi/msep; end end % Compute p-values - pc = (sum(fc<=pdistc)+1)/(arg.nperm+1); - pr = (sum(fr<=pdistr)+1)/(arg.nperm+1); + pc = (sum(fc<=distc)+1)/(arg.nperm+1); + pr = (sum(fr<=distr)+1)/(arg.nperm+1); if reps > 1 - pint = (sum(fi<=pdisti)+1)/(arg.nperm+1); + pint = (sum(fi<=disti)+1)/(arg.nperm+1); p = [pc,pr,pint]; else pint = NaN; @@ -223,12 +218,12 @@ % Compute CIs if nargout > 2 - crit = prctile(pdistc,100*(1-arg.alpha)); + crit = prctile(distc,100*(1-arg.alpha)); cic = [fc./crit;Inf]; - crit = prctile(pdistr,100*(1-arg.alpha)); + crit = prctile(distr,100*(1-arg.alpha)); cir = [fr./crit;Inf]; if reps > 1 - crit = prctile(pdisti,100*(1-arg.alpha)); + crit = prctile(disti,100*(1-arg.alpha)); cii = [fi./crit;Inf]; ci = [cic,cir,cii]; else @@ -268,8 +263,8 @@ % Arrange permutation distributions if nargout > 5 if reps > 1 - pdist = [pdistc,pdistr,pdisti]; + dist = [distc,distr,disti]; else - pdist = [pdistc,pdistr]; + dist = [distc,distr]; end end \ No newline at end of file diff --git a/permutools/permucorr.m b/permutools/permucorr.m index d93da41..36263c4 100644 --- a/permutools/permucorr.m +++ b/permutools/permucorr.m @@ -1,4 +1,4 @@ -function [r,p,ci,stats,pdist] = permucorr(x,varargin) +function [r,p,ci,stats,dist] = permucorr(x,varargin) %PERMUCORR Linear or rank permutation-based correlation. % R = PERMUCORR(X) returns a matrix containing the pairwise linear % correlation coefficients between each pair of columns in X based on @@ -28,7 +28,7 @@ % fields: % 'df' -- the degrees of freedom of each test % -% [R,P,CI,STATS,PDIST] = PERMUCORR(...) returns the permutation +% [R,P,CI,STATS,DIST] = PERMUCORR(...) returns the permuted sampling % distribution of the test statistic. % % [...] = PERMUCORR(...,'PARAM1',VAL1,'PARAM2',VAL2,...) specifies @@ -72,18 +72,20 @@ % PERMUTOOLS https://github.com/mickcrosse/PERMUTOOLS % References: -% [1] Blair RC, Higgins JJ, Karniski W, Kromrey JD (1994) A Study of +% [1] Crosse MJ, Foxe JJ, Molholm S (2024) PERMUTOOLS: A MATLAB +% Package for Multivariate Permutation Testing. arXiv 2401.09401. +% [2] Blair RC, Higgins JJ, Karniski W, Kromrey JD (1994) A Study of % Multivariate Permutation Tests Which May Replace Hotelling's T2 % Test in Prescribed Circumstances. Multivariate Behav Res, % 29(2):141-163. -% [2] Groppe DM, Urbach TP, Kutas M (2011) Mass univariate analysis +% [3] Groppe DM, Urbach TP, Kutas M (2011) Mass univariate analysis % of event-related brain potentials/fields I: A critical tutorial % review. Psychophysiology, 48(12):1711-1725. -% [3] Bishara AJ, Hittner JB, (2012) Testing the Significance of a +% [4] Bishara AJ, Hittner JB, (2012) Testing the Significance of a % Correlation With Nonnormal Data: Comparison of Pearson, % Spearman, Transformation, and Resampling Approaches. Psychol % Methods, 17(3):399-417. -% [4] Bishara AJ, Hittner JB, (2017) Confidence intervals for +% [5] Bishara AJ, Hittner JB, (2017) Confidence intervals for % correlations when data are not normal. Behav Res, 49:294-309. % © 2018-2024 Mick Crosse @@ -173,45 +175,45 @@ [~,idx] = sort(rand(nobs,arg.nperm)); end - % Estimate permutation distribution - pdist = zeros(arg.nperm,nvar); + % Estimate sampling distribution + dist = zeros(arg.nperm,nvar); for i = 1:arg.nperm - pdist(i,:) = (sum(x(idx(:,i),:).*y)-mu)./sdxy; + dist(i,:) = (sum(x(idx(:,i),:).*y)-mu)./sdxy; end % Apply max correction if specified if arg.correct switch arg.tail case 'both' - [~,idx] = max(abs(pdist),[],2); + [~,idx] = max(abs(dist),[],2); csvar = [0;cumsum(ones(arg.nperm-1,1)*nvar)]; - pdist = pdist'; - pdist = pdist(idx+csvar); + dist = dist'; + dist = dist(idx+csvar); case 'right' - pdist = max(pdist,[],2); + dist = max(dist,[],2); case 'left' - pdist = min(pdist,[],2); + dist = min(dist,[],2); end end - % Compute p-value and CIs + % Compute p-value & CI switch arg.tail case 'both' - p = 2*(min(sum(r<=pdist),sum(r>=pdist))+1)/(arg.nperm+1); + p = 2*(min(sum(r<=dist),sum(r>=dist))+1)/(arg.nperm+1); if nargout > 2 - crit = prctile(pdist,100*(1-arg.alpha/2)); + crit = prctile(dist,100*(1-arg.alpha/2)); ci = [max(-1,r-crit);min(1,r+crit)]; end case 'right' - p = (sum(r<=pdist)+1)/(arg.nperm+1); + p = (sum(r<=dist)+1)/(arg.nperm+1); if nargout > 2 - crit = prctile(pdist,100*(1-arg.alpha)); + crit = prctile(dist,100*(1-arg.alpha)); ci = [max(-1,r-crit);Inf(1,nvar)]; end case 'left' - p = (sum(r>=pdist)+1)/(arg.nperm+1); + p = (sum(r>=dist)+1)/(arg.nperm+1); if nargout > 2 - crit = prctile(-pdist,100*(1-arg.alpha)); + crit = prctile(-dist,100*(1-arg.alpha)); ci = [-Inf(1,nvar);min(1,r+crit)]; end end diff --git a/permutools/permuttest.m b/permutools/permuttest.m index 29b219c..389d4f9 100644 --- a/permutools/permuttest.m +++ b/permutools/permuttest.m @@ -1,4 +1,4 @@ -function [t,p,ci,stats,pdist] = permuttest(x,m,varargin) +function [t,p,ci,stats,dist] = permuttest(x,m,varargin) %PERMUTTEST One-sample and paired-sample permutation-based t-test. % T = PERMUTTEST(X) performs a one-sample permutation test based on the % t-statistic of the hypothesis that the data in X come from a @@ -40,7 +40,7 @@ % 'mu' -- the estimated population mean of X, or of X-Y for a % paired test % -% [T,P,CI,STATS,PDIST] = PERMUTTEST(...) returns the permutation +% [T,P,CI,STATS,DIST] = PERMUTTEST(...) returns the permuted sampling % distribution of the test statistic. % % [...] = PERMUTTEST(...,'PARAM1',VAL1,'PARAM2',VAL2,...) specifies @@ -85,13 +85,15 @@ % PERMUTOOLS https://github.com/mickcrosse/PERMUTOOLS % References: -% [1] Blair RC, Higgins JJ, Karniski W, Kromrey JD (1994) A Study of +% [1] Crosse MJ, Foxe JJ, Molholm S (2024) PERMUTOOLS: A MATLAB +% Package for Multivariate Permutation Testing. arXiv 2401.09401. +% [2] Blair RC, Higgins JJ, Karniski W, Kromrey JD (1994) A Study of % Multivariate Permutation Tests Which May Replace Hotelling's T2 % Test in Prescribed Circumstances. Multivariate Behav Res, % 29(2):141-163. -% [2] Gondan M (2010) A permutation test for the race model +% [3] Gondan M (2010) A permutation test for the race model % inequality. Behav Res Methods, 42(1):23-28. -% [3] Groppe DM, Urbach TP, Kutas M (2011) Mass univariate analysis +% [4] Groppe DM, Urbach TP, Kutas M (2011) Mass univariate analysis % of event-related brain potentials/fields I: A critical tutorial % review. Psychophysiology, 48(12):1711-1725. @@ -181,46 +183,46 @@ rng(arg.seed); signx = sign(rand(maxnobs,arg.nperm)-0.5); - % Estimate permutation distribution + % Estimate sampling distribution sqrtn = sqrt(nobs.*df); - pdist = zeros(arg.nperm,nvar); + dist = zeros(arg.nperm,nvar); for i = 1:arg.nperm xp = x.*repmat(signx(:,i),1,nvar); smx = sum(xp,nanflag); - pdist(i,:) = smx./nobs./(sqrt(sum(xp.^2)-(smx.^2)./nobs)./sqrtn); + dist(i,:) = smx./nobs./(sqrt(sum(xp.^2)-(smx.^2)./nobs)./sqrtn); end % Apply max correction if specified if arg.correct - pdist = max(abs(pdist),[],2); + dist = max(abs(dist),[],2); end % Add negative values - pdist(arg.nperm+1:2*arg.nperm,:) = -pdist; + dist(arg.nperm+1:2*arg.nperm,:) = -dist; arg.nperm = 2*arg.nperm; if arg.verbose fprintf('Adding negative of values to permutation distribution.\n') fprintf('Number of permutations used: %d\n',arg.nperm) end - % Compute p-value and CIs + % Compute p-value & CI switch arg.tail case 'both' - p = 2*(sum(abs(t)<=pdist)+1)/(arg.nperm+1); + p = 2*(sum(abs(t)<=dist)+1)/(arg.nperm+1); if nargout > 2 - crit = prctile(pdist,100*(1-arg.alpha/2)).*se; + crit = prctile(dist,100*(1-arg.alpha/2)).*se; ci = [mu-crit;mu+crit]; end case 'right' - p = (sum(t<=pdist)+1)/(arg.nperm+1); + p = (sum(t<=dist)+1)/(arg.nperm+1); if nargout > 2 - crit = prctile(pdist,100*(1-arg.alpha)).*se; + crit = prctile(dist,100*(1-arg.alpha)).*se; ci = [mu-crit;Inf(1,nvar)]; end case 'left' - p = (sum(t>=pdist)+1)/(arg.nperm+1); + p = (sum(t>=dist)+1)/(arg.nperm+1); if nargout > 2 - crit = prctile(pdist,100*(1-arg.alpha)).*se; + crit = prctile(dist,100*(1-arg.alpha)).*se; ci = [-Inf(1,nvar);mu+crit]; end end diff --git a/permutools/permuttest2.m b/permutools/permuttest2.m index a2ef9fb..741a4b7 100644 --- a/permutools/permuttest2.m +++ b/permutools/permuttest2.m @@ -1,4 +1,4 @@ -function [t,p,ci,stats,pdist] = permuttest2(x,y,varargin) +function [t,p,ci,stats,dist] = permuttest2(x,y,varargin) %PERMUTTEST2 Unpaired two-sample permutation-based t-test. % T = PERMUTTEST2(X,Y) performs a two-sample permutation test based on % the t-statistic of the hypothesis that the data in X and Y come from @@ -37,7 +37,7 @@ % deviations (unequal variance) % 'mu' -- the estimated population mean of X-Y % -% [T,P,CI,STATS,PDIST] = PERMUTTEST2(...) returns the permutation +% [T,P,CI,STATS,DIST] = PERMUTTEST2(...) returns the permuted sampling % distribution of the test statistic. % % [...] = PERMUTTEST2(...,'PARAM1',VAL1,'PARAM2',VAL2,...) specifies @@ -76,17 +76,19 @@ % PERMUTOOLS https://github.com/mickcrosse/PERMUTOOLS % References: -% [1] Blair RC, Higgins JJ, Karniski W, Kromrey JD (1994) A Study of +% [1] Crosse MJ, Foxe JJ, Molholm S (2024) PERMUTOOLS: A MATLAB +% Package for Multivariate Permutation Testing. arXiv 2401.09401. +% [2] Blair RC, Higgins JJ, Karniski W, Kromrey JD (1994) A Study of % Multivariate Permutation Tests Which May Replace Hotelling's T2 % Test in Prescribed Circumstances. Multivariate Behav Res, % 29(2):141-163. -% [2] Groppe DM, Urbach TP, Kutas M (2011) Mass univariate analysis +% [3] Groppe DM, Urbach TP, Kutas M (2011) Mass univariate analysis % of event-related brain potentials/fields I: A critical tutorial % review. Psychophysiology, 48(12):1711-1725. -% [3] Groppe DM, Urbach TP, Kutas M (2011) Mass univariate analysis +% [4] Groppe DM, Urbach TP, Kutas M (2011) Mass univariate analysis % of event-related brain potentials/fields II: Simulation % studies. Psychophysiology, 48(12):1726-1737. -% [4] Groppe DM (2016) Combating the scientific decline effect with +% [5] Groppe DM (2016) Combating the scientific decline effect with % confidence (intervals). Psychophysiology, 54(1):139-145. % © 2018-2024 Mick Crosse @@ -185,8 +187,8 @@ i1 = idx(1:maxnobsx,:); i2 = idx(maxnobsx+1:maxnobs,:); - % Estimate permutation distribution - pdist = zeros(arg.nperm,nvar); + % Estimate sampling distribution + dist = zeros(arg.nperm,nvar); for i = 1:arg.nperm x1 = x(i1(:,i),:); x2 = x(i2(:,i),:); @@ -200,36 +202,36 @@ case 'unequal' sep = sqrt(var1./nobsx+var2./nobsy); end - pdist(i,:) = (sm1./nobsx-sm2./nobsy)./sep; + dist(i,:) = (sm1./nobsx-sm2./nobsy)./sep; end % Apply max correction if specified if arg.correct - [~,idx] = max(abs(pdist),[],2); + [~,idx] = max(abs(dist),[],2); csvar = [0;cumsum(ones(arg.nperm-1,1)*nvar)]; - pdist = pdist'; - pdist = pdist(idx+csvar); + dist = dist'; + dist = dist(idx+csvar); end - % Compute p-value and CIs + % Compute p-value & CI switch arg.tail case 'both' - pdabs = abs(pdist); + pdabs = abs(dist); p = (sum(abs(t)<=pdabs)+1)/(arg.nperm+1); if nargout > 2 crit = prctile(pdabs,100*(1-arg.alpha)).*se; ci = [mu-crit;mu+crit]; end case 'right' - p = (sum(t<=pdist)+1)/(arg.nperm+1); + p = (sum(t<=dist)+1)/(arg.nperm+1); if nargout > 2 - crit = prctile(pdist,100*(1-arg.alpha)).*se; + crit = prctile(dist,100*(1-arg.alpha)).*se; ci = [mu-crit;Inf(1,nvar)]; end case 'left' - p = (sum(t>=pdist)+1)/(arg.nperm+1); + p = (sum(t>=dist)+1)/(arg.nperm+1); if nargout > 2 - crit = prctile(pdist,100*(1-arg.alpha)).*se; + crit = prctile(dist,100*(1-arg.alpha)).*se; ci = [-Inf(1,nvar);mu+crit]; end end diff --git a/permutools/permuvartest2.m b/permutools/permuvartest2.m index 9332369..f215d26 100644 --- a/permutools/permuvartest2.m +++ b/permutools/permuvartest2.m @@ -1,4 +1,4 @@ -function [f,p,ci,stats,pdist] = permuvartest2(x,y,varargin) +function [f,p,ci,stats,dist] = permuvartest2(x,y,varargin) %PERMUVARTEST2 Unpaired two-sample permutation-based F-test. % F = PERMUVARTEST2(X,Y) performs a two-sample permutation test based on % the F-statistic of the null hypothesis that the data in X and Y come @@ -27,7 +27,7 @@ % 'df1' -- the numerator degrees of freedom of each test % 'df2' -- the denominator degrees of freedom of each test % -% [F,P,CI,STATS,PDIST] = PERMUVARTEST2(...) returns the permutation +% [F,P,CI,STATS,DIST] = PERMUVARTEST2(...) returns the permuted sampling % distribution of the test statistic. % % [...] = PERMUVARTEST2(...,'PARAM1',VAL1,'PARAM2',VAL2,...) specifies @@ -62,11 +62,13 @@ % PERMUTOOLS https://github.com/mickcrosse/PERMUTOOLS % References: -% [1] Blair RC, Higgins JJ, Karniski W, Kromrey JD (1994) A Study of +% [1] Crosse MJ, Foxe JJ, Molholm S (2024) PERMUTOOLS: A MATLAB +% Package for Multivariate Permutation Testing. arXiv 2401.09401. +% [2] Blair RC, Higgins JJ, Karniski W, Kromrey JD (1994) A Study of % Multivariate Permutation Tests Which May Replace Hotelling's T2 % Test in Prescribed Circumstances. Multivariate Behav Res, % 29(2):141-163. -% [2] Groppe DM, Urbach TP, Kutas M (2011) Mass univariate analysis +% [3] Groppe DM, Urbach TP, Kutas M (2011) Mass univariate analysis % of event-related brain potentials/fields I: A critical tutorial % review. Psychophysiology, 48(12):1711-1725. @@ -145,32 +147,32 @@ i1 = idx(1:maxnobsx,:); i2 = idx(maxnobsx+1:maxnobs,:); - % Estimate permutation distribution - pdist = zeros(arg.nperm,nvar); + % Estimate sampling distribution + dist = zeros(arg.nperm,nvar); for i = 1:arg.nperm x1 = x(i1(:,i),:); x2 = x(i2(:,i),:); var1 = (sum(x1.^2,nanflag)-(sum(x1,nanflag).^2)./nobsx)./df1; var2 = (sum(x2.^2,nanflag)-(sum(x2,nanflag).^2)./nobsy)./df2; - pdist(i,:) = var1./var2; + dist(i,:) = var1./var2; end % Apply max correction if specified if arg.correct - [~,imax] = max(pdist,[],2); - [~,imin] = min(pdist,[],2); + [~,imax] = max(dist,[],2); + [~,imin] = min(dist,[],2); csvar = [0;cumsum(ones(arg.nperm-1,1)*nvar)]; - pdist = pdist'; - pdmax = pdist(imax+csvar); - pdmin = pdist(imin+csvar); + dist = dist'; + pdmax = dist(imax+csvar); + pdmin = dist(imin+csvar); k = 1; else - pdmax = pdist; - pdmin = pdist; + pdmax = dist; + pdmin = dist; k = 2; end - % Compute p-value and CIs + % Compute p-value & CI switch arg.tail case 'both' p = k*(min(sum(f<=pdmax),sum(f>=pdmin))+1)/(arg.nperm+1); diff --git a/permutools/permuztest.m b/permutools/permuztest.m index 80ae8c6..24b8359 100644 --- a/permutools/permuztest.m +++ b/permutools/permuztest.m @@ -1,4 +1,4 @@ -function [z,p,ci,stats,pdist] = permuztest(x,m,sigma,varargin) +function [z,p,ci,stats,dist] = permuztest(x,m,sigma,varargin) %PERMUZTEST One-sample permutation-based Z-test. % Z = PERMUZTEST(X,M,SIGMA) performs a one-sample permutation test based % on the Z-statistic of the null hypothesis that the data in X come from @@ -25,7 +25,7 @@ % 'sd' -- the estimated population standard deviation of X % 'mu' -- the estimated population mean of X % -% [Z,P,CI,STATS,PDIST] = PERMUZTEST(...) returns the permutation +% [Z,P,CI,STATS,DIST] = PERMUZTEST(...) returns the permuted sampling % distribution of the test statistic. % % [...] = PERMUZTEST(...,'PARAM1',VAL1,'PARAM2',VAL2,...) specifies @@ -63,11 +63,13 @@ % PERMUTOOLS https://github.com/mickcrosse/PERMUTOOLS % References: -% [1] Blair RC, Higgins JJ, Karniski W, Kromrey JD (1994) A Study of +% [1] Crosse MJ, Foxe JJ, Molholm S (2024) PERMUTOOLS: A MATLAB +% Package for Multivariate Permutation Testing. arXiv 2401.09401. +% [2] Blair RC, Higgins JJ, Karniski W, Kromrey JD (1994) A Study of % Multivariate Permutation Tests Which May Replace Hotelling's T2 % Test in Prescribed Circumstances. Multivariate Behav Res, % 29(2):141-163. -% [2] Groppe DM, Urbach TP, Kutas M (2011) Mass univariate analysis +% [3] Groppe DM, Urbach TP, Kutas M (2011) Mass univariate analysis % of event-related brain potentials/fields I: A critical tutorial % review. Psychophysiology, 48(12):1711-1725. @@ -118,47 +120,47 @@ rng(arg.seed); signx = sign(rand(maxnobs,arg.nperm)-0.5); - % Estimate permutation distribution + % Estimate sampling distribution diffxm = x-m; sen = se.*nobs; - pdist = zeros(arg.nperm,nvar); + dist = zeros(arg.nperm,nvar); for i = 1:arg.nperm xp = diffxm.*repmat(signx(:,i),1,nvar); smx = sum(xp,nanflag); - pdist(i,:) = smx./sen; + dist(i,:) = smx./sen; end % Apply max correction if specified if arg.correct - pdist = max(abs(pdist),[],2); + dist = max(abs(dist),[],2); end % Add negative values - pdist(arg.nperm+1:2*arg.nperm,:) = -pdist; + dist(arg.nperm+1:2*arg.nperm,:) = -dist; arg.nperm = 2*arg.nperm; if arg.verbose fprintf('Adding negative of values to permutation distribution.\n') fprintf('Number of permutations used: %d\n',arg.nperm) end - % Compute p-value and CIs + % Compute p-value & CI switch arg.tail case 'both' - p = 2*(sum(abs(z)<=pdist)+1)/(arg.nperm+1); + p = 2*(sum(abs(z)<=dist)+1)/(arg.nperm+1); if nargout > 2 - crit = prctile(pdist,100*(1-arg.alpha/2)).*se; + crit = prctile(dist,100*(1-arg.alpha/2)).*se; ci = [mu-crit;mu+crit]; end case 'right' - p = (sum(z<=pdist)+1)/(arg.nperm+1); + p = (sum(z<=dist)+1)/(arg.nperm+1); if nargout > 2 - crit = prctile(pdist,100*(1-arg.alpha)).*se; + crit = prctile(dist,100*(1-arg.alpha)).*se; ci = [mu-crit;Inf(1,nvar)]; end case 'left' - p = (sum(z>=pdist)+1)/(arg.nperm+1); + p = (sum(z>=dist)+1)/(arg.nperm+1); if nargout > 2 - crit = prctile(pdist,100*(1-arg.alpha)).*se; + crit = prctile(dist,100*(1-arg.alpha)).*se; ci = [-Inf(1,nvar);mu+crit]; end end