Efficient and safe Monte Carlo testing with multiple-comparisons correction

MMCTest - Monte-carlo Multiple Comparison testing:
--- Overview

This function performs simultaneous multiple Monte Carlo tests to estimate significance over a set of hypotheses, while controlling either the false discovery rate or the familywise error rate. The algorithm is an implementation of Gandy & Hahn 2014 [1]; notation closely follows this publication.

--- Basic usage
Usage: [vbRejectH0] = MMCTest(fhMCSample, nm)
[vbRejectH0, vbAcceptH0, vbUndecided, mfCI, vnk] = ...
MMCTest(fhMCSample, nm <, fAlpha, nc, nk_max, fEpsilon, fhCI, fhNeta, fhMCHTest>)

'fhMCSample' is a Monte Carlo sampling function (either with or without replacement), with the signature
mbExceedTestStat = @fhMCSample(vnHypotheses, nNumSamples)

This function must return an [M x N] matrix, where each element is a single sample from one hypothesis 'vnHypothesis(m)' under test. Each element is a boolean variable that indicates whether the MC resampled test statistic for the corresponding hypothesis was equal to or more extreme that the observed test statistic for that hypothesis. I.e., 'mbExceedTestStat' is equivalent to X_ij in [1]. Sampling can be performed with or without replacement, as desired, and is controlled entirely wihtin 'fhMCSample'. Elements of 'vnHypotheses' are [1..nm].

'nm' is a scalar integer specifying the number of hypotheses under test.

'vbRejectH0' will be a [1 x nm] boolean vector, indicating whether H0 should be rejected for the corresponding hypothesis. By default, the test will be performed at a threshold of 'fAlpha' = 0.05.

The optional argument 'fAlpha' can be used to specify the test significance threshold.


--- Return arguments

'vbRejectH0', 'vbAcceptH0' and 'vbUndecided' are [1 x nm] boolean vectors, specifying whether H0 should be rejected, accepted, or could not be decided after a fixed number of samples at the 'fAlpha' significance level.

'mfCI' is a [2 x nm] matrix, containing confidence intervals for P for each hypothesis. Each [2 x 1] column vector is ['fUpperPCI'; 'fLowerPCI'].

'vnk' is a [1 x nm] integer vector, containing the number of samples taken from each corresponding hypothesis.


--- Example usage

% - Construct a hypothesis testing function, with a true P value of 0.06
>> fhMCSample = @(vnH, nNS)rand(numel(vnH), nNS) < 0.06;

% - Perform a test with two hypotheses, with default parameters (i.e. fAlpha = 0.05).
>> [vbRejectH0, vbAcceptH0, vbUndecided, mfCI, vnk] = MMCTest(fhMCSample, 2)

vbRejectH0 =
0 0

vbAcceptH0 =
1 1

vbUndecided =
0 0

mfCI =
0.0520 0.0488
0.0801 0.0769

vnk =
5566 5566

H0 is correctly accepted for both hypotheses at fAlpha = 0.05. 'vnk' indicates that 5566 Monte Carlo samples were drawn for each hypothesis. Both confidence intervals ('mfCI') contain the true P value of 0.06.

% - Perform a test with fAlpha = 0.061
>> [vbRejectH0, ~, ~, mfCI, vnk] = MMCTest(fhMCSample, 2, 0.061)
vbRejectH0 =
1 1

mfCI =
0.0594 0.0592
0.0610 0.0606

vnk =
2872983 2872983

H0 is correctly rejected for both hypotheses at fAlpha = 0.061. Many more Monte Carlo samples were required to reach this decision, since 'fAlpha' was close to the true P value. Both confidence intervals contain the true P value of 0.06.

--- References

[1] Gandry & Hahn 2014, "MMCTest - A Safe Algorithm for Implementing Multiple Monte Carlo Tests". Scandinavian J. of Statistics 41, pp1083--1101. DOI: 10.1111/sjos.12085