123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117 |
- % Bonferroni-Holm (1979) correction for multiple comparisons. This is a
- % sequentially rejective version of the simple Bonferroni correction for multiple
- % comparisons and strongly controls the family-wise error rate at level alpha.
- %
- % It works as follows:
- % 1) All p-values are sorted in order of smallest to largest. m is the
- % number p-values.
- % 2) If the 1st p-value is greater than or equal to alpha/m, the procedure
- % is stopped and no p-values are significant. Otherwise, go on.
- % 3) The 1st p-value is declared significant and now the second p-value is
- % compared to alpha/(m-1). If the 2nd p-value is greater than or equal
- % to alpha/(m-1), the procedure is stopped and no further p-values are
- % significant. Otherwise, go on.
- % 4) Et cetera.
- %
- % As stated by Holm (1979) "Except in trivial non-interesting cases the
- % sequentially rejective Bonferroni test has strictly larger probability of
- % rejecting false hypotheses and thus it ought to replace the classical
- % Bonferroni test at all instants where the latter usually is applied."
- %
- %
- % function [corrected_p, h]=bonf_holm(pvalues,alpha)
- %
- % Required Inputs:
- % pvalues - A vector or matrix of p-values. If pvalues is a matrix, it can
- % be of any dimensionality (e.g., 2D, 3D, etc...).
- %
- % Optional Input:
- % alpha - The desired family-wise alpha level (i.e., the probability of
- % rejecting one of more null hypotheses when all null hypotheses are
- % really true). {default: 0.05}
- %
- % Output:
- % corrected_p - Bonferroni-Holm adjusted p-values. Any adjusted p-values
- % less than alpha are significant (i.e., that null hypothesis
- % is rejected). The adjusted value of the smallest p-value
- % is p*m. The ith smallest adjusted p-value is the max of
- % p(i)*(m-i+1) or adj_p(i-1). Note, corrected p-values can
- % be greater than 1.
- % h - A binary vector or matrix of the same dimensionality as
- % pvalues. If the ith element of h is 1, then the ith p-value
- % of pvalues is significant. If the ith element of h is 0, then
- % the ith p-value of pvalues is NOT significant.
- %
- % Example:
- % >>p=[.56 .22 .001 .04 .01]; %five hypothetical p-values
- % >>[cor_p, h]=bonf_holm(p,.05)
- % cor_p =
- % 0.5600 0.4400 0.0050 0.1200 0.0400
- % h =
- % 0 0 1 0 1
- %
- % Conclusion: the third and fifth p-values are significant, but not the
- % remaining three.
- %
- % Reference:
- % Holm, S. (1979) A simple sequentially rejective multiple test procedure.
- % Scandinavian Journal of Statistics. 6, 65-70.
- %
- %
- % For a review on contemporary techniques for correcting for multiple
- % comparisons that are often more powerful than Bonferroni-Holm see:
- %
- % Groppe, D.M., Urbach, T.P., & Kutas, M. (2011) Mass univariate analysis
- % of event-related brain potentials/fields I: A critical tutorial review.
- % Psychophysiology, 48(12) pp. 1711-1725, DOI: 10.1111/j.1469-8986.2011.01273.x
- % http://www.cogsci.ucsd.edu/~dgroppe/PUBLICATIONS/mass_uni_preprint1.pdf
- %
- %
- % Author:
- % David M. Groppe
- % Kutaslab
- % Dept. of Cognitive Science
- % University of California, San Diego
- % March 24, 2010
- function [corrected_p, h]=bonf_holm(pvalues,alpha)
- if nargin<1,
- error('You need to provide a vector or matrix of p-values.');
- else
- if ~isempty(find(pvalues<0,1)),
- error('Some p-values are less than 0.');
- elseif ~isempty(find(pvalues>1,1)),
- fprintf('WARNING: Some uncorrected p-values are greater than 1.\n');
- end
- end
- if nargin<2,
- alpha=.05;
- elseif alpha<=0,
- error('Alpha must be greater than 0.');
- elseif alpha>=1,
- error('Alpha must be less than 1.');
- end
- s=size(pvalues);
- if isvector(pvalues),
- if size(pvalues,1)>1,
- pvalues=pvalues';
- end
- [sorted_p sort_ids]=sort(pvalues);
- else
- [sorted_p sort_ids]=sort(reshape(pvalues,1,prod(s)));
- end
- [dummy, unsort_ids]=sort(sort_ids); %indices to return sorted_p to pvalues order
- m=length(sorted_p); %number of tests
- mult_fac=m:-1:1;
- cor_p_sorted=sorted_p.*mult_fac;
- cor_p_sorted(2:m)=max([cor_p_sorted(1:m-1); cor_p_sorted(2:m)]); %Bonferroni-Holm adjusted p-value
- corrected_p=cor_p_sorted(unsort_ids);
- corrected_p=reshape(corrected_p,s);
- h=corrected_p<alpha;
|