PERMUTOOLS/examples/run_booteffectsize_examples.m at master · mickcrosse/PERMUTOOLS · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
function run_booteffectsize_examples
%RUN_BOOTEFFECTSIZE_EXAMPLES  Run bootstrapped effect size examples.
%   Generates random multivariate data for 2 samples X and Y. Each sample
%   has 20 variables, each with a mean value of 0, except for the first 10
%   variables of Y which have a mean value of -1. Each variable has 30
%   observations. Effect sizes based on Cohen's d and Glass' Δ with
%   bootstrapped confidence intervals (CIs) are measured between the
%   corresponding variables of each sample for independent and dependent
%   samples, as well as samples with equal and unequal variances. The
%   results are compared to those of the equivalent parametric measures
%   (i.e. CIs using the Student's t-distribution) using meanEffectSize.m.
%
%   Note: The parametric effect sizes and CIs generated by meanEffectSize.m
%   for the uncorrected plots are actually bias-corrected (as there is no
%   option to specify bias-corrected in meanEffectSize.m) and are thus
%   slightly overestimated.
%
%   See also BOOTEFFECTSIZE MEANEFFECTSIZE.
%
%   PERMUTOOLS https://github.com/mickcrosse/PERMUTOOLS

%   References:
%       [1] Crosse MJ, Foxe JJ, Molholm S (2024) PERMUTOOLS: A MATLAB
%           Package for Multivariate Permutation Testing. arXiv 2401.09401.

%   © 2018-2026 Mick Crosse <crossemj@tcd.ie>
%   CNL, Albert Einstein College of Medicine, NY.
%   TCBE, Trinity College Dublin, Ireland.

% Set up experiment
nobs = 30; nvar = 20;
xaxis = 1:nvar;
paired = [false,true];
samples = {'indep.','dep.'};
vartype = {'equal','unequal'};

% Generate random data
rng(42);
x = randn(nobs,nvar);
y = randn(nobs,nvar);
y(:,1:round(nvar/2)) = y(:,1:round(nvar/2))-1;

% Set up figure
figure('Name','Effect size analysis: effect sizes & CIs','NumberTitle','off')
set(gcf,'color','w')

k = 1;
for i = 1:numel(paired)

    for n = 1:numel(vartype)

        % Make variance unequal
        if n == 2
            yp = y*1.25;
        else
            yp = y;
        end

        % Set effect size measure
        if paired(i)==true && strcmp(vartype{n},'unequal')
            paired(i) = false;
            samples{i} = 'indep.';
            effect = 'Glass';
        else
            effect = 'Cohen';
        end

        % Plot parametric & uncorrected bootstrapped CIs
        d1 = zeros(1,nvar);
        ci1 = zeros(2,nvar);
        for j = 1:nvar
            stats1 = meanEffectSize(x(:,j),yp(:,j),'Effect',effect,...
                'Paired',paired(i),'VarianceType',vartype{n},...
                'ConfidenceIntervalType','exact');
            d1(j) = stats1.Effect;
            ci1(:,j) = stats1.ConfidenceIntervals';
        end
        [d2,ci2] = booteffectsize(x,yp,'paired',paired(i),...
            'vartype',vartype{n},'effect',effect,'correct',0);
        subplot(4,2,k), hold on
        plot(xaxis,d2,'LineWidth',3)
        plot(xaxis,ci1,'k',xaxis,ci2,'--r')
        xlim([0,nvar+1]), ylim([-2,5]), box on, grid on
        if i == 1 && n == 1
            title('Uncorrected')
        end
        if i == 2 && n == 2
            xlabel('variable')
        end
        ylabel({[samples{i},' samples'];['w/ ',vartype{n},' SD']})
        switch effect
            case 'Cohen'
                legend('Cohen''s {\itd}','95% CI (param.)','',...
                    '95% CI (boot.)','Location','best')
            case 'Glass'
                legend('Glass'' {\itΔ}','95% CI (param.)','',...
                    '95% CI (boot.)','Location','best')
        end
        [d2,ci2] = booteffectsize(x,yp,'paired',paired(i),...
            'vartype',vartype{n},'effect',effect,'correct',1);
        subplot(4,2,k+1), hold on
        plot(xaxis,d2,'LineWidth',3)
        plot(xaxis,ci1,'k',xaxis,ci2,'--r')
        xlim([0,nvar+1]), ylim([-2,5]), box on, grid on
        if i == 1 && n == 1
            title('Bias-corrected')
        end
        if i == 2 && n == 2
            xlabel('variable')
        end
        switch effect
            case 'Cohen'
                legend('Hedges'' {\itg}','95% CI (param.)','',...
                    '95% CI (boot.)','Location','best')
            case 'Glass'
                legend('Glass'' {\itΔ}','95% CI (param.)','',...
                    '95% CI (boot.)','Location','best')
        end

        k = k+2;

    end

end