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.

To compute, however, the exact p-value which would produce specific p-values of form p= rather than p<alpha (for assessing sensitivity to see how close a p-value's unadjusted version is to statistical significance) this [attachment:fdrexact.xls spreadsheet] can be used.

Reference

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.

[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

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.

Lix LM and Sajobi T (2010) Testing multiple outcomes in repeated measures designs. Psychological Methods 15(3) 268-280.