Skip to content

Commit

Permalink
Renamed dist var to avoid overwrite
Browse files Browse the repository at this point in the history
  • Loading branch information
@mickcrosse
mickcrosse committed Jan 18, 2024
1 parent 1f90341 commit e2ae758
Show file tree
Hide file tree
Showing 8 changed files with 142 additions and 158 deletions.
44 changes: 14 additions & 30 deletions permutools/booteffectsize.m
View file
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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 <[email protected]>
Expand Down Expand Up @@ -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]));
Expand All @@ -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);
Expand All @@ -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);
Expand All @@ -348,8 +339,6 @@
diffxyb(:,j) = diffi(:);
end
end

% Compute mean difference
switch arg.effect
case 'cliff'
mub = sum(diffxyb,nanflag)/nobs(i);
Expand All @@ -358,8 +347,6 @@
otherwise
mub = smxb/nobsx(i)-smyb/nobsy(i);
end

% Compute standard deviation
switch arg.effect
case 'glass'
sdb = sqrt(varxb);
Expand All @@ -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

Expand Down
27 changes: 11 additions & 16 deletions permutools/permuanova1.m
View file
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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 <[email protected]>
% CNL, Albert Einstein College of Medicine, NY.
Expand Down Expand Up @@ -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);
Expand All @@ -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

Expand Down
45 changes: 20 additions & 25 deletions permutools/permuanova2.m
View file
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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 <[email protected]>
% CNL, Albert Einstein College of Medicine, NY.
Expand Down Expand Up @@ -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));
Expand All @@ -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;
Expand All @@ -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
Expand Down Expand Up @@ -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
44 changes: 23 additions & 21 deletions permutools/permucorr.m
View file
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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 <[email protected]>
Expand Down Expand Up @@ -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
Expand Down
Loading

0 comments on commit e2ae758

Please sign in to comment.