SPSS script to Perform the Benjamini and Hochberg FDR procedure
This procedure may also be computed using the [:FAQ/pvs:spreadsheet here]. Keselman HJ, Cribbie R and Holland B (1999) compared several pairwise comparison tests and found the FDR approach to be the most powerful when more than 6 pairwise tests were performed e.g. when comparing all possible pairs of means for five groups. Benjamini and Hochberg (1995) also show FDR has higher power compared to other procedures including a p-value adjustment procedure due to Hochberg (1988) which, in turn, is shown as having higher power than both Bonferroni and Holm methods for multiple testing in repeated measures by Lix and Sajobi (2010).
The FDR (Benjamini and Hochberg, 1995) is a stepdown procedure using a prescribed alpha usually 0.05 for testing in advance of the p-value testing and so the above spreadsheet should be used for hypothesis testing of multiple p-values. See [:FAQ/SpssBonferroni: here] for a warning about using exact p-values instead of statements of form p<alpha.
The exact p-value which produces specific p-values of form p= rather than p<alpha may, however, be computed using this [attachment:fdrexact.xls spreadsheet]. The exact p-value could be computed to assess sensitivity to see how close an adjusted p-value is to statistical significance.
References
Benjamini Y and Hochberg Y (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society B 57 289-300.
Hochberg Y (1988) A sharper Bonferroni procedure for multiple tests of significance. Biometrika 75 800-803.
Keselman HJ, Cribbie R and Holland B (1999) The pairwise multiple comparison multiplicity problem: an alternative approach to familywise and comparisonwise type I error control Psychological methods 4(1) 58-69.
Lix LM and Sajobi T (2010) Testing multiple outcomes in repeated measures designs. Psychological Methods 15(3) 268-280.
[COPY AND PASTE THE BELOW BOX SYNTAX INTO A SPSS SYNTAX WINDOW; SELECT ALL AND RUN. EDIT THE INPUT DATA AS REQUIRED]
DATA LIST FREE / pvals .
BEGIN DATA
0.021
0.001
0.017
0.041
0.005
0.036
0.042
0.023
0.07
0.1
END DATA .
COMPUTE alpha = 0.05 .
COMPUTE index = $CASENUM .
EXECUTE .
SORT CASES BY pvals (A) .
COMPUTE rank = index .
COMPUTE ref = alpha * rank / 10 .
COMPUTE diff = pvals - ref .
COMPUTE comp = diff LT 0 .
CREATE ccomp = CSUM(comp) .
EXECUTE .
MATRIX .
GET rnk / VARIABLES = rank .
GET pvals / VARIABLES = pvals .
GET ccomp / VARIABLES = ccomp .
GET alpha / VARIABLES = alpha .
GET diff / VARIABLES = diff .
COMPUTE ccmpmx = pvals LE CMAX(((CMAX(ccomp) &* (diff LE 0)) EQ ccomp) &* pvals) .
SAVE {rnk, pvals, alpha, ccmpmx} / OUTFILE = * / VARIABLES = index,pvals,alpha,fdrsig .
END MATRIX .
SORT CASES BY index (A) .
EXPORT OUTFILE = out .
IMPORT FILE = out / KEEP = pvals alpha fdrsig .
FORMAT pvals (F5.3) .
FORMAT alpha (F5.3) .
FORMAT fdrsig (F1.0) .
VALUE LABELS fdrsig 1 'Sig' 0 'Nonsig' .
SET RESULTS = LISTING .
DO IF $CASENUM EQ 1 .
PRINT EJECT /'P Value' 1 'FDR Criterion' 9 'FDR Significance' 24 .
END IF .
PRINT / pvals (T3 F5.3) alpha (T17, F5.3) fdrsig (T39, F1.0) .
EXECUTE .This produces the following output:
P Value FDR Criterion FDR Significance .021 .050 1 .001 .050 1 .017 .050 1 .041 .050 0 .005 .050 1 .036 .050 0 .042 .050 0 .023 .050 1 .070 .050 0 .100 .050 0
R Code to Perform the Benjamini and Hochberg FDR procedure
[Insert P values in any order] pvalue <- c(0.021,0.001,0.017,0.041,0.005,0.036,0.042,0.023,0.07,0.1) sorted.pvalue<-sort(pvalue) sorted.pvalue j.alpha <- (1:10)*(.05/10) diff <- sorted.pvalue-j.alpha neg.diff <- diff[diff<0] pos.diff <- neg.diff[length(neg.diff)] index <- diff==pos.diff p.cutoff <-sorted.pvalue[index] p.cutoff p.sig <- pvalue[pvalue <= p.cutoff] p.sig
R Output
> pvalue<-c(0.021,0.001,0.017,0.041,0.005,0.036,0.042,0.023,0.07,0.1) > sorted.pvalue <- sort(pvalue) > sorted.pvalue [sorted P values in ascending order] [1] 0.001 0.005 0.017 0.021 0.023 0.036 0.041 0.042 0.070 0.100 > j.alpha <- (1:10) * (0.05/10) > diff <- sorted.pvalue - j.alpha > neg.diff <- diff[diff < 0] > pos.diff <- neg.diff[length(neg.diff)] > index <- diff == pos.diff > p.cutoff <- sorted.pvalue[index] > p.cutoff [1] 0.023 > p.sig <- pvalue[pvalue <= p.cutoff] > p.sig [significant p values by B-H method]: [1] 0.021 0.001 0.017 0.005 0.023